CN102609946A - 一种基于黎曼流形的脑白质纤维束跟踪的组间处理方法 - Google Patents

Abstract

本发明公开了一种基于黎曼流形的脑白质纤维束跟踪的组间处理方法，该方法包括以下步骤：将脑白质看作黎曼流形，并分别在Fick第二定律和黎曼流形下讨论扩散及各向异性，再通过Navier-Stoke方程表达扩散张量场中流体的运动，从而构造一个偏微分方程(PDE)来描述黎曼流形下两点间测地线的性质，该偏微分方程的数值解通过水平集(Level-Set)方法计算得到。在纤维束跟踪过程中采用一种感兴趣区(ROI)逆配准方法，消除了个体之间脑解剖结构的差异性，最后进行纤维束数目的组间处理。本发明利用整个张量来构造黎曼流形，充分利用了张量场的所有信息，并通过PDE和Level-Set方法很好地解决了抗噪性问题，提高了鲁棒性。

Description

一种基于黎曼流形的脑白质纤维束跟踪的组间处理方法

技术领域

本发明属于磁共振扩散张量成像(Diffusion Tensor Imaging，DTI)领域，具体涉及一种基于黎曼(Riemannian)流形的脑白质纤维束跟踪，再对纤维束的数目进行组间处理的方法。

背景技术

扩散张量成像(DTI)是一种磁共振技术，是目前唯一可以无创地对脑组织中水分子的各向异性扩散性质进行定量分析和测量的技术。DTI通过测量活体组织中水分子毫米数量级的扩散运动速率与状态，反映水分子的各向异性，从而得到纤维组织的走向，表现人体组织的微观结构。近年来，扩散张量成像的研究主要集中在脑成像方面，特别是应用于脑白质纤维束的成像。黎曼流形是一个微分流形，其中每点p的切空间都定义了点积，而且其数值随p平滑地改变。每个Rn的平滑子流形可以导出黎曼度量：把Rn的点积都限制于切空间内。我们可以定义黎曼流形为和Rn的平滑子流形是等距同构的度量空间，等距是指其内蕴度量(intrinsic metric)和上述从Rn导出的度量是相同的。虽然黎曼流形通常是弯曲的，但“直线”的概念依然存在：那就是测地线。

近年来，DTI研究的重要进展是将其应用于脑白质纤维的可视化，其中脑白质纤维束跟踪是DTI可视化研究的一个热点。在脑白质纤维中，沿着长轴方向扩散的水分子由于阻力较小而扩散速度较快，而垂直于长轴方向扩散的水分子由于阻力较大而扩散速度较慢。纤维跟踪使用连续的曲线表示纤维的走向和分布，通过纤维跟踪方法可以得到三维连续的脑白质组织结构，并可以显示脑白质纤维细节。传统的脑白质纤维束成像算法大多将脑白质中的各个体素用一个3×3的张量来表示，并计算各个张量的最大特征值和最大特征向量，以最大特征向量的方向来表示该体素的纤维束主要走向，从而形成一个三维的矢量场。

目前对脑白质神经纤维的跟踪方法已经有很多种，大概可以分成基于张量域和基于全局能量最小化的两种方法。基于张量域的纤维跟踪算法主要利用局部张量信息进行纤维跟踪，DTI能产生每个体素的优选扩散方向，空间上每个张量的排列称为张量域。在该三维矢量场中选取一个种子点，再运用插值、波传导等方法迭代计算源自该种子点的脑白质纤维束，迭代算法的终止条件一般为纤维束传导至各向异性较低的区域，或纤维束相邻体素间偏折角度过大等。基于张量域的纤维跟踪算法关键在于当前点扩散方向的确定，但其存在一个缺点，不能处理纤维分叉和交叉。

而基于全局能量最小化的方法可以处理分叉和交叉纤维，该方法增加了纤维不确定性和随机性的考虑，通过最小化代价函数，寻求最优路径，具有最小代价的路径与真实路径相对应，但处理过程中引入了太多的假阳性支，并且对于纤维走向描述不明确。

因此，传统的脑白质纤维跟踪方法具有抗噪性差、纤维束跟踪的鲁棒性低的缺点。并且，由于不同人的脑解剖结构存在显著差异，传统的脑白质纤维跟踪方法在比较不同人的纤维束特征时会有困难。

发明内容

(一)要解决的技术问题

为了克服已有技术的不足，本发明提供一种基于黎曼流形的脑白质纤维束跟踪的组间处理方法。

(二)技术方案

本发明提出一种基于黎曼流形脑白质纤维束跟踪的组间处理方法，包括以下步骤：

步骤Sa：对脑白质DTI图像数据进行配准校正，在配准校正后的脑白质DTI图像数据上建立Fick第二定理公式；

步骤Sb：为脑白质DTI图像数据构造黎曼流形，选取两点间测地线作为脑白质纤维束的路径，其中两点间测地线即为最小距离路径；

步骤Sc：计算所述黎曼流形下的扩散算子；

步骤Sd：依据所述扩散算子，构造Navier-Stoke方程表达扩散张量场中流体的运动，从而构造一个偏微分方程来描述黎曼流形下两点间测地线的性质；

步骤Se：依据Level-Set数值解方法，通过迭代计算所述偏微分方程，演化得出脑白质纤维束成像结果；

步骤Sf：对所述脑白质纤维束成像结果使用感兴趣区逆配准的方法确定脑白质纤维束跟踪的初始区域，得到特定的纤维束；

步骤Sg，使用个体样本的所述特定脑白质纤维束的数目进行处理分析，得到组间处理结果。

本发明采用基于黎曼流形的脑白质纤维束成像算法进行脑白质纤维束的跟踪。脑白质在扩散张量成像形成的张量场中，每个体素都可由一个对称、正定、协变的张量表示，在此基础上将脑白质看作黎曼流形，并分别在Fick第二定律和黎曼流形下讨论扩散及各向异性，再通过Navier-Stoke方程表达扩散张量场中流体的运动，从而构造一个偏微分方程(PDE)来描述黎曼流形下两点间测地线的性质，该偏微分方程的数值解通过Level-Set方法计算得到。

(三)有益效果

本发明利用整个张量来构造黎曼流形，充分利用了张量场的所有信息，并通过PDE和Level-Set方法很好地解决了抗噪性问题，提高了鲁棒性。在纤维束跟踪过程中采用一种感兴趣区(ROI)逆配准方法跟踪不同分组的脑白质纤维束，消除了个体之间脑解剖结构的差异性，可进行纤维束数目处理上的组间处理，为不同个体或不同组之间的脑白质纤维束的准确比较提供一种可行的方法。同时，本发明为脑白质微结构的分析和可视化提供了一种有效途径，不仅有助于深入地了解脑纤维的结构，还可为纤维缺失或结构异常造成的疾病诊断提供有效的信息。

附图说明

图1为本发明基于黎曼流形的脑白质纤维束跟踪的组间处理方法的流程示意图。

图2是本发明实施例中使用所述基于脑白质纤维束跟踪的组间处理方法获得的各向异性分数(FA)图及有颜色编码的FA图。

图3是本发明实施例中使用所述基于脑白质纤维束跟踪的组间处理方法获得的彩色FA图，黑色箭头所指的地方(蓝色)为内囊后肢。

图4是本发明实施例中使用所述基于脑白质纤维束跟踪的组间处理方法获得的皮质脊髓束(CST)脑白质纤维图，A行表示一个正常人，B行表示一个中风病人。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

参照图1，本发明涉及一种基于脑白质纤维束跟踪的组间处理方法，尤其涉及一种基于黎曼流形和Level-Set方法由扩散张量成像(DTI)跟踪出脑白质纤维束，再对纤维束的数量进行组间处理的方法，其具体实施步骤如下：

步骤Sa，对脑白质DTI图像数据进行配准校正，为配准校正的脑白质DTI图像建立Fick第二定理公式；

1.对脑白质DTI图像进行配准校正

由于磁共振扫描过程中各种各样的噪声的影响，个体自身存在尺度和位置上的差异，非常有必要在分析数据之前对数据做一定的预处理。在整个的实验的数据获取中，主要的噪声信息来源有：(1)物理头动；(2)图像内层间扫描时间差别；(3)外在磁场的不均匀性等。DTI预处理是在保留脑功能图像细节的同时，使用非共线加权方向的图像与非加权像(B₀像)进行仿射配准变换预处理，并去掉一些质量不好方向的图像，提高脑功能图像的信噪比。非共线加权方向的图像与非加权像是使用传统的自旋回波-弥散张量成像技术(SE-DTI)获取的。这里的仿射配准变换是指对各个非共线加权方向的图像进行缩放、旋转、平移后，根据非共线加权方向的图像和非加权像间的共有特征进行几何配准。

2.为配准校正的DTI图像数据建立Fick第二定理公式

布朗运动是一种无规则的运动，如水中的水分子在一直不停地做着布朗运动。但在生物体中，由于水分子在不同微环境下具有不同的渗透能力，水分子的布朗运动可能被独立的大分子、细胞内器官等抑制。生物体中的这种现象即为各向异性。这种扩散的各向异性可以通过扩散加权成像(Diffusion Weighted Imaging，DWI)来测量。如果在图像数据采集时得到6个或6个以上不同方向的DWI图像数据，就可以得到图像数据中任意一点的扩散张量D，并以3×3对称矩阵表示：

$D=\left[\begin{array}{ccc}{D}_{\mathrm{xx}}& D_{\mathrm{xy}}& D_{\mathrm{xz}}\\ D_{\mathrm{xy}}& D_{\mathrm{yy}}& D_{\mathrm{yz}}\\ D_{\mathrm{xz}}& D_{\mathrm{yz}}& D_{\mathrm{zz}}\end{array}\right]$

其中x，y，z表示三维空间坐标的三个方向，D表示扩散张量，D_xx，D_xy，D_xz表示沿着x轴、沿xoy平面和xoz平面方向的扩散权值。

D可以通过解以下线性方程组得到：

$\left(\begin{array}{cccccc}& & & \·& & \\ & & & \·& & \\ & & & \·& & \\ {x_i}^2& {y_i}^2& {z_i}^2& 2x_i{y}_i& 2y_i{z}_i& 2z_i{x}_i\\ & & & \·& & \\ & & & \·& & \\ & & & \·& & \end{array}\right)\left(\begin{array}{c}{D}_{\mathrm{xx}}\\ D_{\mathrm{yy}}\\ D_{\mathrm{zz}}\\ D_{\mathrm{xy}}\\ D_{\mathrm{yz}}\\ D_{\mathrm{zx}}\end{array}\right)=\left(\begin{array}{c}\·\\ \·\\ \·\\ -\frac{1}{b}\ln \frac{S_i}{S_0}\\ \·\\ \·\\ \·\end{array}\right)$

其中i＝1，…，n，n表示扫描方向数(n≥6)，(x_i，y_i，z_i)^T表示扩散梯度脉冲方向；S_i表示施加了扩散梯度脉冲得到的扩散加权图像数据；S₀表示未施加扩散梯度脉冲的图像数据；b表示扩散敏感系数。

菲克(Fick)定律是描述扩散规律性的定律。由德国生理学家菲克(1829-1901)于1855年提出。菲克(Fick)定律包括两个内容：

(1)扩散第一定律。***中i物质的扩散达到稳定状态时，也即i物质在各处的浓度分别不随时间而变化时，单位时间内通过垂直于扩散方向单位横截面的i物质的通量J_i与i物质的浓度梯度

(x方向上单位距离上i物质的浓度差，其中C_i表示i物质的浓度)成正比，即

式中D_i为比例常数，亦称“扩散系数”，其量纲为L²·T^-1(L为长度，T为时间)。因

总是负值，故在等式右方加负号，以使J_i为正值。

(2)扩散第二定律。扩散过程尚未达到稳定状态前，i物质浓度C_i随时间t和位置(只考虑x方向)而变化的关系，服从偏微分式

对于具体的扩散过程，要利用其特定的起始条件和边界条件求解此式，得出C_i＝f(x，t)的具体函数。该定律是处理各种扩散传质过程理论的有力工具。

由Fick第二定律可知，随着时间的增加，扩散会对浓度场造成影响，由Fick第二定理导出下式：

$u_t=\▿\·(D\▿u)$

其中表示梯度，t表示时间，D表示DTI图像数据的扩散张量，u表示扩散浓度。

步骤Sb，为脑白质DTI图像数据构造黎曼流形，选取两点间测地线(最小距离路径)作为脑白质纤维束的路径；

扩散张量成像中的扩散张量D是对称、正定、协变的，因此，可将G＝D^-1(其中上标-1表示求逆运算)作为黎曼度量，从而将脑白质转化为一个黎曼流形，并将

所表示的扩散算子与黎曼度量相联系，使计算黎曼流形下的测地线路径、距离与脑白质纤维束相联系。测地线是黎曼流形上两点之间的局部最短距离，是黎曼流形的内蕴几何特征。测地线距离定义如下：设(M，G)是一个黎曼流形，x、y分别为M上的两点，则x到y的测地线距离为M上连接x、y两点的最短测地线长度，记为d(x，y)，

$d(x,y)={\arg \min }_{\mathrm{\γ}\∈\mathrm{\Γx},y}{\∫}_0^{T_{\mathrm{xy}}}{\left|{\mathrm{\γ}}^{\′}\left(t\right)\right|}_R\mathrm{dt}$

其中，Γx，y表示连接x、y两点的测地线，γ(0)＝x，γ(T_xy)＝y，且|γ′(t)|_R＝1，arg表示目标函数的参数，min表示求最小值，T_xy表示使得γ函数值为y的一个常量，t表示时间参数，dt表示对时间求导。

通过G＝D^-1(其中上标-1表示求逆运算)可知，DTI中特征值较大的张量所对应的测地线距离较小，因此，传统脑白质纤维束成像算法中选取最大特征向量作为纤维束主要走向，对应在黎曼流形中即为选取两点间测地线(最小距离路径)作为纤维束的路径。

步骤Sc，计算所述黎曼流形下的扩散算子；

黎曼流形(M，G)下扩散算子为

${\▿}_G^{2}u=G^{\mathrm{ij}}(\frac{{\∂}^{2}u}{\∂x^k\∂x^l}-{\mathrm{\Γ}}^l\frac{\∂u}{\∂x^l})=\frac{1}{\sqrt{\left|G\right|}}\frac{\∂}{\∂x^k}\left(\sqrt{\left|G\right|}{G}^{\mathrm{kl}}\frac{\∂u}{\∂x^l}\right)$

其中G＝D^-1为扩散张量矩阵D的逆矩阵，u表示扩散浓度，x表示某一坐标方向，

表示求二阶梯度，

表示在u表达式中对x^l求偏导数，

表示在u表达式中先对x^k求偏导数，再对x^l求偏导数，

表示在u表达式中对x^k求偏导数，G^ij，G^kl表示扩散张量矩阵D的逆矩阵的第i行j列元素和第k行l列元素，其中1≤i，j，k，l≤3，且均为整数。将此扩散算子与步骤Sa中的Fick第二定理公式联立得

$\left\{\begin{array}{c}{u}_t=\▿\·(D\▿u)\\ {\▿}_G^{2}u=\frac{1}{\sqrt{\left|G\right|}}\frac{\∂}{\∂x^k}\left(\sqrt{\left|G\right|}{G}^{\mathrm{kl}}\frac{\∂u}{\∂x^l}\right)\end{array}\right.$

解以上方程组可得

$\&dtri;\&CenterDot;(D\&dtri;u)={\&dtri;}_G^{2}u-\frac{1}{2}<\&dtri;\lg \left|G\right|,{\&dtri;}_{G}u>$

其中

表示梯度，lg表示取10为底的对数，<.>表示求内积。

步骤Sd，依据上一步骤确定的扩散算子，构造Navier-Stoke方程表达扩散张量场中流体的运动，从而构造一个偏微分方程(PDE)来描述所述黎曼流形下两点间测地线的性质；

Navier-Stoke方程是描述液体和空气这样的流体物质运动的方程，其建立了流体的粒子动量加速度和作用在液体内部的压力的改变和粘滞力以及引力之间的关系，这些粘滞力产生于分子的相互作用，使得该方程可以描述作用于液体任意给定区域的力的动态平衡。Navier-Stoke方程可以模拟扩散及脑白质纤维束走向：

$\mathrm{\&rho;}(\frac{\&PartialD;v}{\&PartialD;t}+v\&CenterDot;\&dtri;v)=\&dtri;\&CenterDot;G-\&dtri;p$

其中，v表示流速；ρ表示流体密度；G表示扩散张量场；p表示压力，

表示在v表达式中对t求偏导数，

表示梯度。

同时，在光滑、连通、完备的黎曼流形(M，G)中，测地线长度函数u为标量函数，令K为M的一个闭包、非空子集，则有：

$|\&dtri;u|=1,u\left(x\right)=0,x\&Element;K,K\&Subset;M$

将此关系式代入Navier-Stoke方程，并联立步骤Sc中的扩散算子方程得：

$\left\{\begin{array}{c}\mathrm{\&rho;}\frac{\&PartialD;v}{\&PartialD;t}={\&dtri;}_G^{2}u-\frac{1}{2}<\&dtri;\lg \left|G\right|,{\&dtri;}_{G}u>-\&dtri;p-\mathrm{\&rho;v}\&CenterDot;\&dtri;v\\ v_0=0\\ \frac{\&PartialD;u}{\&PartialD;N}{|}_{\&PartialD;\mathrm{\&Omega;}}=0,\&ForAll;u\&Element;\&PartialD;\mathrm{\&Omega;}\end{array}\right.$

其中

表示求二阶梯度，

表示梯度，lg表示取10为底的对数，<.>表示求内积，

表示在v表达式中对t求偏导数，

表示在u表达式中对N求偏导数，Ω表示待处理的扩散张量图像数据，

表示脑白质边界条件，

表示任意。上述方程组的第1个等式可以写成

p(v，v′)+Q(u，u′)＝0

上式即为一个偏微分方程(PDE)，其中p(.)，Q(.)表示两个函数，v′为v的微分，u′为u的微分。可以通过哈密顿-雅可比(Hamilton-Jacobi)方程的数值解法来求解。Hamilton-Jacobi方程的一般形式如下，

$\left\{\begin{array}{c}\frac{\&PartialD;v}{\&PartialD;t}+H(x_1,\&CenterDot;\&CenterDot;\&CenterDot;,x_n,t,\frac{\&PartialD;v}{\&PartialD;x_1},\&CenterDot;\&CenterDot;\&CenterDot;,\frac{\&PartialD;v}{\&PartialD;x_n})=0,\\ v(x_1,\&CenterDot;\&CenterDot;\&CenterDot;,x_n,0)=v_0(x_1,\&CenterDot;\&CenterDot;\&CenterDot;,x_n).\end{array}\right.$

其中

表示在v表达式中对t求偏导数，H(.)表示Hamilton函数，

(i＝1，...，n)表示在v表达式中对x_i求偏导数，n表示空间维数，v₀(x₁，…，x_n)表示v表达式的初值。这类方程来源于最优控制理论、微分几何、物理学、流体力学等。Hamilton-Jacobi方程解的性质非常复杂。一方面，H一J方程的弱解存在但不唯一；另一方面，具有连续的解但其导数却不连续。

步骤Se，依据Level-Set数值解方法，通过迭代计算所述偏微分方程(PDE)演化得出脑白质纤维束成像结果；

Level-Set数值解方法可解决应用于图像数据或运动分割的曲线或曲面演化间题，其基本数学思想是：将当前正在演化的曲线看作更高一维函数的水平集，利用曲线或曲面演化与Hamilton-Jacobi方程的相似性，给出一种准确、鲁棒的计算方法。Level-Set函数在演化过程中始终保持为一个函数，因此，可以在一个时空离散的网格中表示Level-Set函数。设网格间隔为1，每个节点用(i，j，k)表示，离散化的时间用t_n＝nr，n∈N表示，r为时间步长，节点(i，j，k)在时间t_n的Level-Set函数值可以用

来表示。通过引入有限差分算子及有限差分，函数P、Q可以分别在时空离散网格中表示为：

$P=\mathrm{\&rho;}[\frac{v_{i,j,k}^{n+1}-v_{i,j,k}^n}{r}+\max (v_{i,j,k}^n,0){\&dtri;}^{+}+\min (v_{i,j,k}^n,0){\&dtri;}^{-}]$

$Q=g^{11}[\max {(D_{i,j,k}^{x-}u,0)}^2+\max {(D_{i,j,k}^{x+}u,0)}^2]+g^{22}[\max {(D_{i,j,k}^{y-}u,0)}^2+\max {(D_{i,j,k}^{y+}u,0)}^2]$

$+g^{33}[\max {(D_{i,j,k}^{z-}u,0)}^2+\max {(D_{i,j,k}^{z+}u,0)}^2]+2g^{12}[\max (D_{i,j,k}^{x-}u,0)+\min (D_{i,j,k}^{x+}u,0)]\&CenterDot;$

$[\max (D_{i,j,k}^{y-}u,0)+\max (D_{i,j,k}^{y+}u,0)]+2g^{13}[\max (D_{i,j,k}^{x-}u,0)+\min (D_{i,j,k}^{x+}u,0)]\&CenterDot;$

$[\max (D_{i,j,k}^{z-}u,0)+\max (D_{i,j,k}^{z+}u,0)]+2g^{23}[\max (D_{i,j,k}^{y-}u,0)+\min (D_{i,j,k}^{y+}u,0)]\&CenterDot;$

$[\max (D_{i,j,k}^{z-}u,0)+\max (D_{i,j,k}^{z+}u,0)]$

其中g^ij∈G＝D^-1，为扩散张量矩阵的元素，max表示求最大值，min表示求最小值。由此通过迭代计算演化得出脑白质纤维束成像结果。

步骤Sf，对所述脑白质纤维束成像结果使用ROI逆配准的方法确定脑白质纤维束跟踪的初始区域，得到特定的脑白质纤维束；

特定的纤维束例如是如皮质脊髓束(corticospinal tract，CST)纤维束，感兴趣区(ROI)逆配准方法是在标准脑模板上做一个标准ROI，与每个被试的非加权像进行配准，得到每个被试的ROI模板。

根据国际标准大脑模板EPI做一个ROI的解剖结构保存为mask，即所选择的ROI区域体素坐标置为1，周围其余体素坐标全是0。利用仿射变换配准方法把从标准模板生成的ROI解剖结构的mask配准到每个被试的非加权B0像上，形成新的相同解剖结构的mask，这样克服了个体间的解剖结构差异。所谓仿射配准变换方式是指进行对原始图像缩放、旋转、平移后，根据原始图像和标准模板间的共有特征进行几何配准。

有各向异性分数(fractional anisotropy，FA)是用来定量分析各向异性的参数，是水分子各向异性成分占整个弥散张量的比例，它的变化范围从0～1。0代表弥散不受限制，比如脑脊液的FA值接近0；对于非常规则的具有方向性的组织，其FA值大于0，例如脑白质纤维FA值接近1。我们设定一个FA阈值，如0.2，在所选ROI内大于等于此阈值的所有体素被作为跟踪的起始种子区域，利用前述的步骤Sb-Se的方法跟踪出纤维束，并处理出每组个体纤维束的数目。

步骤Sg，使用个体样本的所述特定脑白质纤维束的数目进行组间处理分析，得到组间处理结果，即两组整体的特定脑白质纤维束数目的均值差，作为比较两组样本的特定脑白质纤维束差异的一个指标。

使用样本数据生成经验假设分布的方法进行处理检验。对两组纤维束数目进行双样本t检验，得到一个p′值。T检验，亦称student t检验(Student′s t test)，主要用于样本含量较小(例如n＜30)，总体标准差未知的正态分布资料。T检验是用于小样本(样本容量小于30)的两个平均值差异程度的检验方法。p′值是将观察结果认为有效即具有总体代表性的犯错概率，在这里即使第一组的纤维束数目的均值与第二组的纤维束数目的均值不相等，也会在p′值的概率下相等，即允许有一定误差的意思，犯错的概率为p′。

本发明所述的基于黎曼流形的脑白质纤维束跟踪的组间处理方法，可通过真实的扩散张量成像数据得以说明：

(1)真实数据实验

为展示本发明的效果，在实施方案中采用真实数据集作测试，共12个被测试者参与了实验，分为两组，一组是左侧大脑损伤的偏瘫病人，4男2女，另一组是正常人，也是4男2女。实验过程由专业的影像科机师操作，实验中采用单次触发自旋回波平面回波成像(single-shotspin-echo echo-planar imaging)序列获取DTI图像数据。

采用SPM5(http://www.fil.ion.ucl.ac.uk/spm/)进行标准EPI模板的ROI解剖结构(左右内囊后肢)选取及对非加权的B0像ROI数据逆配准。然后将脑白质看作黎曼流形，并分别在Fick第二定律和黎曼流形下讨论扩散及各向异性，再通过Navier-Stoke方程表达扩散张量场中流体的运动，从而构造一个偏微分方程(PDE)来描述黎曼流形下两点间测地线的性质，该偏微分方程的数值解通过Level-Set方法计算得到。得到全脑的纤维束跟踪结果之后，将之前所选的ROI作为跟踪起始区域对每个被测试者的皮质脊髓束(corticospinal tract，CST)纤维束数目进行组间处理。

(2)实验结果

在真实实验数据集上采用本发明的方法进行纤维束跟踪，并用逆配准的方法进行CST数目组间处理比较。如图2和图3所示为基于黎曼流形分解出的各向异性分数(FA)图及有颜色编码的FA图，并由此图可准确确定左右大脑的内囊后肢为ROI。用在标准EPI模板上确定的内囊后肢mask对正常人和中风病人进行逆配准，得到消除个体解剖差异的ROI，作为个体纤维束跟踪的初始区域。图4为获得的皮质脊髓束(CST)白质纤维图。A行表示一个正常人，B行表示一个中风病人。从图上就可以观察出中风病人比正常人的纤维束有明显的缺失。我们对每组的纤维束数量进行两样本t检验，得到表格1与表格2，两组左侧纤维束数目的均值具有显著差异(p′＝0.0001)，而右侧无明显差异(p′＝0.8230)。

以上实验结果说明，本发明所述的一种基于黎曼流形的脑白质纤维束跟踪，再对纤维束的数量进行组间处理的方法，可以有效地消除个体之间脑解剖结构的差异性，进行纤维束数目处理上的组间处理，为不同组之间的脑白质纤维束的准确比较提供一种可行的方法。

表格1正常组与中风组的左侧CST纤维束数目双样本检验

表格2正常组与中风组的右侧CST纤维束数目双样本检验

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (9)

1.一种基于黎曼流形的脑白质纤维束跟踪的组间处理方法，其特征在于，包括以下步骤：

步骤Sa：对脑白质DTI图像数据进行配准校正，在配准校正后的脑白质DTI图像数据上建立Fick第二定理公式；

步骤Sb：为脑白质DTI图像数据构造黎曼流形，选取两点间测地线作为脑白质纤维束的路径，其中两点间测地线即为最小距离路径；

步骤Sc：计算所述黎曼流形下的扩散算子；

步骤Sd：依据所述扩散算子，构造Navier-Stoke方程表达扩散张量场中流体的运动，从而构造一个偏微分方程来描述黎曼流形下两点间测地线的性质；

步骤Se：依据Level-Set数值解方法，通过迭代计算所述偏微分方程，演化得出脑白质纤维束成像结果；

步骤Sf：对所述脑白质纤维束成像结果使用感兴趣区逆配准的方法确定脑白质纤维束跟踪的初始区域，得到特定的纤维束；

步骤Sg，使用个体样本的所述特定脑白质纤维束的数目进行处理分析，得到组间处理结果。

2.如权利要求1所述的基于黎曼流形的脑白质纤维束跟踪的组间处理方法，其特征在于，所述DTI图像数据预处理是在保留DTI图像数据细节的同时，使非共线加权方向的图像数据与非加权B0像进行仿射配准变换预处理。

3.如权利要求1所述的基于黎曼流形的脑白质纤维束跟踪的组间处理方法，其特征在于，利用Fick第二定律可知，对扩散张量图像数据进行时间上的描述，如下式所示：

$u_t=\&dtri;\&CenterDot;(D\&dtri;u)$

其中

表示梯度，t表示时间，D表示扩散张量，u表示扩散浓度。

4.如权利要求3所述的基于黎曼流形的脑白质纤维束跟踪的组间处理方法，其特征在于构造黎曼流形，选取两点间测地线(最小距离路径)作为纤维束的路径，其步骤为；

A、定义测地线距离：设(M，G)是一个黎曼流形，x、y分别为M上的两点，则x到y的测地线距离为M上连接x、y两点的最短测地线长度，记为d(x，y)，

$d(x,y)={\arg \min }_{\mathrm{\&gamma;}\&Element;\mathrm{\&Gamma;x},y}{\&Integral;}_0^{T_{\mathrm{xy}}}{\left|{\mathrm{\&gamma;}}^{\&prime;}\left(t\right)\right|}_R\mathrm{dt}$

其中，Γx，y表示连接x、y两点的测地线，γ(0)＝x，γ(T_xy)＝y，且|γ′(t)|_R＝1，arg表示目标函数的参数，min表示求最小值。

B、通过G＝D^-1(其中上标-1表示求逆运算)可知，DTI图像中特征值较大的张量所对应的测地线距离较小，因此，传统脑白质纤维束成像算法中选取最大特征向量作为纤维束主要走向，对应在黎曼流形中即为选取两点间测地线(最小距离路径)作为纤维束的路径。

5.如权利要求1所述的基于黎曼流形的脑白质纤维束跟踪的组间处理方法，其特征在于，计算黎曼流形下的扩散算子；

A、黎曼流形(M，G)下扩散算子为

${\&dtri;}_G^{2}u=G^{\mathrm{ij}}(\frac{{\&PartialD;}^{2}u}{\&PartialD;x^k\&PartialD;x^l}-{\mathrm{\&Gamma;}}^l\frac{\&PartialD;u}{\&PartialD;x^l})=\frac{1}{\sqrt{\left|G\right|}}\frac{\&PartialD;}{\&PartialD;x^k}\left(\sqrt{\left|G\right|}{G}^{\mathrm{kl}}\frac{\&PartialD;u}{\&PartialD;x^l}\right)$

其中G＝D^-1为扩散张量矩阵D的逆矩阵，u表示扩散浓度，x表示某一坐标方向，

表示求二阶梯度，

表示在u表达式中对x^l求偏导数，

表示在u表达式中先对x^k求偏导数，再对x^l求偏导数。

B、将此扩散算子与步骤Sa中的Fick第二定理公式联立得

$\left\{\begin{array}{c}{u}_t=\&dtri;\&CenterDot;(D\&dtri;u)\\ {\&dtri;}_G^{2}u=\frac{1}{\sqrt{\left|G\right|}}\frac{\&PartialD;}{\&PartialD;x^k}\left(\sqrt{\left|G\right|}{G}^{\mathrm{kl}}\frac{\&PartialD;u}{\&PartialD;x^l}\right)\end{array}\right.$

解以上方程组可得

$\&dtri;\&CenterDot;(D\&dtri;u)={\&dtri;}_G^{2}u-\frac{1}{2}<\&dtri;\lg \left|G\right|,{\&dtri;}_{G}u>$

其中表示梯度，lg表示取10为底的对数，<.>表示求内积。

6.如权利要求1所述的基于黎曼流形的脑白质纤维束跟踪的组间处理方法，其特征在于，步骤Sd具体包括以下步骤；

A、Navier-Stoke方程可以模拟扩散及脑白质纤维束走向：

$\mathrm{\&rho;}(\frac{\&PartialD;v}{\&PartialD;t}+v\&CenterDot;\&dtri;v)=\&dtri;\&CenterDot;G-\&dtri;p$

其中，v表示流速；ρ表示流体密度；G表示扩散张量场；p表示压力，表示在v表达式中对t求偏导数，表示梯度。

B、在光滑、连通、完备的黎曼流形(M，G)中，测地线长度函数u为标量函数，令K为M的一个闭包、非空子集，则有：

$|\&dtri;u|=1,u\left(x\right)=0,x\&Element;K,K\&Subset;M$

C、将此关系式代入Navier-Stoke方程，并联立步骤Sc中的扩散算子方程得：

$\left\{\begin{array}{c}\mathrm{\&rho;}\frac{\&PartialD;v}{\&PartialD;t}={\&dtri;}_G^{2}u-\frac{1}{2}<\&dtri;\lg \left|G\right|,{\&dtri;}_{G}u>-\&dtri;p-\mathrm{\&rho;v}\&CenterDot;\&dtri;v\\ v_0=0\\ \frac{\&PartialD;u}{\&PartialD;N}{|}_{\&PartialD;\mathrm{\&Omega;}}=0,\&ForAll;u\&Element;\&PartialD;\mathrm{\&Omega;}\end{array}\right.$

其中

表示求二阶梯度，表示梯度，lg表示取10为底的对数，<.>表示求内积，

表示在v表达式中对t求偏导数，

表示在u表达式中对N求偏导数，Ω表示脑白质边界条件，

表示任意。

D、上述方程组的第1个等式可以写成

p(v，v′)+Q(u，u′)＝0

其中p(.)，Q(.)表示两个函数，v′为v的微分，u′为u的微分。

7.如权利要求6所述的基于黎曼流形的脑白质纤维束跟踪的组间处理方法，其特征在于，依据Level-Set数值解方法，通过迭代计算演化得出脑白质纤维束成像结果，其步骤为；

A、Level-Set函数在演化过程中始终保持为一个函数，因此，可以在一个时空离散的网格中表示Level-Set函数。设网格间隔为1，每个节点用(i，j，k)表示，离散化的时间用表示t_n＝nr，n∈N，r为时间步长，节点(i，j，k)在时间t_n的Level-Set函数值可以用

来表示。

B、通过引入有限差分算子及迎风有限差分，函数P、Q可以分别在时空离散网格中表示为：

$P=\mathrm{\&rho;}[\frac{v_{i,j,k}^{n+1}-v_{i,j,k}^n}{r}+\max (v_{i,j,k}^n,0){\&dtri;}^{+}+\min (v_{i,j,k}^n,0){\&dtri;}^{-}]$

$Q=g^{11}[\max {(D_{i,j,k}^{x-}u,0)}^2+\max {(D_{i,j,k}^{x+}u,0)}^2]+g^{22}[\max {(D_{i,j,k}^{y-}u,0)}^2+\max {(D_{i,j,k}^{y+}u,0)}^2]$

$+g^{33}[\max {(D_{i,j,k}^{z-}u,0)}^2+\max {(D_{i,j,k}^{z+}u,0)}^2]+2g^{12}[\max (D_{i,j,k}^{x-}u,0)+\min (D_{i,j,k}^{x+}u,0)]\&CenterDot;$

$[\max (D_{i,j,k}^{y-}u,0)+\max (D_{i,j,k}^{y+}u,0)]+2g^{13}[\max (D_{i,j,k}^{x-}u,0)+\min (D_{i,j,k}^{x+}u,0)]\&CenterDot;$

$[\max (D_{i,j,k}^{z-}u,0)+\max (D_{i,j,k}^{z+}u,0)]+2g^{23}[\max (D_{i,j,k}^{y-}u,0)+\min (D_{i,j,k}^{y+}u,0)]\&CenterDot;$

$[\max (D_{i,j,k}^{z-}u,0)+\max (D_{i,j,k}^{z+}u,0)]$

其中g^ij∈G＝D^-1，为扩散张量矩阵的元素，max表示求最大值，min表示求最小值。由此通过迭代计算演化得出脑白质纤维束成像结果。

8.如权利要求7所述的基于黎曼流形的脑白质纤维束跟踪的组间处理方法，其特征在于，步骤Sf具体的步骤为；

A、根据国际标准大脑模板EPI做一个ROI的解剖结构保存为mask。

B、利用仿射变换配准方法把从标准模板生成的ROI解剖结构的mask配准到每个被试的非加权B0像上，形成新的相同解剖结构的mask。

C、设定一个FA阈值，在所选ROI内大于等于此阈值的所有体素被作为跟踪的起始种子区域，利用前述的步骤Sb-Se的方法跟踪出的纤维束。

D、处理出每组个体纤维束的数目。

9.如权利要求1所述的基于黎曼流形的脑白质纤维束跟踪的组间处理方法，其特征在于，步骤Sg组分析使用样本数据生成经验假设分布的方法进行处理检验。对两组纤维束数目进行双样本t检验，得到一个p′值。可比较两组数据DTI脑白质纤维束数目的差异是否显著。

Images