一种非负高阶张量拟牛顿搜索的纤维方向分布估计方法
【技术领域】
[0001] 本发明涉及医学成像、图像处理、数值分析、三维重建、计算机科学、神经解剖学等 领域，尤其是一种用于脑白质纤维成像的非负纤维方向分布的高阶笛卡尔扩散张量成像方 法。
【背景技术】
[0002] 脑白质纤维成像是一种非侵入式的获取人脑白质区域的体素神经纤维微结构信 息并进行三维重构和展示的信息技术。其主要方法是对原始的扩散加权磁共振数据进行体 素建模，获得每个体素内的纤维方向分布情况，形成具有解剖学意义的纤维空间微结构。通 过大量研究，该领域的学者在脑纤维微结构重构算法、不确定信息处理等方面已取得了一 系列成果。然而，当前流行的纤维重构算法大都无法保证所求得的纤维方向分布函数的非 负性。这显然不符合扩散信号的实际物理意义，同时，它还会在一定程度上促成伪峰的形 成，从而影响纤维成像的结果，极大地阻碍了该技术在临床医学中的应用。
【发明内容】
[0003] 为了克服现有纤维成像方法无法保证扩散函数非负性的问题，本发明提出一种以 高阶张量为导向的基于拟牛顿搜索的非负纤维方向分布估计方法。
[0004] 本发明所采用的技术方案如下所示：
[0005] -种非负高阶张量拟牛顿搜索的纤维方向分布估计方法，所述的纤维方向分布估 计方法包括以下步骤：
[0006] (1)数据预处理:读取脑部扩散加权磁共振数据，获取施加梯度方向g时的磁共振 信号S(g)和未施加梯度方向时的磁共振信号So,以及相应的梯度方向数据，选取感兴趣区 域，并计算该区域的扩散衰减信号S(g)/S〇;
[0007] (2)将感兴趣区域内每个体素中的扩散衰减信号逐个建模为具有扩散形态的椭球 分布模型，建模过程如下：
[0008] 2.1)体素微结构建模:将扩散衰减信号假设为沿重建向量v的单条纤维信号响应 函数R (v，g)与扩散函数D (v)在球面上的卷积：
[0009]
[0010] 其中，糾,.4) = 1?-近似看作一个高斯分布函数，g= {giERlx31 i = l，· · ·，n}为梯 度方向，v={vPeRlx3|p = l，. . .，K}为在单位球面上采样的重建向量，R1X3表示维度为1X3 的实数域矩阵空间，η和K分别表示梯度方向和重建向量的个数，μ = ε b是表征扩散效率ε与 扩散敏感系数b共同影响的一个参数，扩散函数D(v)的表达式如下：
[0011]
[0012 ] drs表示单项式v;V:. V丨ts的系数，1为高阶张量的阶数，r，s分别表示重建向量v = (Vx， Vy，VZ)的基方向VX和Vy的指数，表示第j个张量的张量系数，j = l，. . .，m，m表示张量的个数
，./? = ν:? -表示第j个张量单项式，并满足p + + + 1，c是由m个张 量系数组成的系数向量，表达式为〇=|>1，\2，.人]1^(>) = [；^(>)彳2(>)，...；^(>)]由1]1个 张量单项式构成；
[0013] 2.2)数学模型：
[0014] 扩散加权磁共振信号有η个扩散梯度方向gi，i = l，...，n，并且沿重建向量v进行 重建，那么系数向量c通过最小化下面的代价函数J(c)求得：
[0015]
[0016] 其中，Ei = S(gi)/So是第i个扩散梯度方向gi上的衰减信号； 是一个mXm维的矩阵，其值只与扩散梯度方向gl、重建向量v以及参数μ有关，对于每一个扩 散梯度方向gi，都有一个Qi矩阵与之对应；
[0017] (3)计算张量系数向量c，得到扩散函数D(v)，再计算每个采样点处的扩散函数值， 最后将扩散函数值拟合成扩散模型，搜索极值并计算纤维方向。
[0018] 进一步，所述步骤(3)中，所述张量系数向量c的计算包括以下步骤：
[0019] 3.1)在单位半球面上均匀采样321个离散的点，以球心为原点获取这321个重建向 量v，计算单条纤维响应函数R(v，g)的值，设定高阶张量模型的阶数1，计算单项式矩阵F (v)，进而计算出步骤2.2)中的矩阵&;
[0020] 3.2)使用BFGS拟牛顿搜索算法迭代求解2.2)中的最小化问题，步骤如下：
[0021]步骤3.2.1已知代价函数J(c)，选取一个初始点(^作为第一次迭代的搜索起始点， 迭代次数计为k=l，并设置最大迭代次数k_max，计算J(c)的梯度向量表达式：
[0022]
[0023] 步骤3.2.2对于第k次迭代，令表示第k次迭代的梯度向量，计算拟牛顿方
向dk = _Hktk作为第k次迭代的搜索方向，以ck为起点，沿方向dk进行一维搜索，求得本次搜索 的可接受步长ak，矩阵Hk表示的是第k次迭代时⑶点处Hesse矩阵的逆的近似矩阵，其更新方 法如下-
[0024]
[0025] 其中，I是单位阵，Sk = Ck+1-Ck是相邻两次迭代的解的差，0k = tk+1_tk是相邻两次迭 代的梯度向量之差；
[0026] 步骤3.2.3更新张量系数向量Ck+1 = Ck+akdk，计算对应的代价函数值J(ck+1)，当达 到最大迭代次数或相邻两次迭代的代价函数值满足时，则终止迭代，否则， 返回步骤3.2.2，进入下一次迭代;参数σ为一个较小的常数，用于判断算法是否收敛到一个 局部极小值。
[0027]再进一步，所述步骤(3)还包括3.3)将得到的张量系数用于拟合扩散函数，获取纤 维方向分布函数模型，搜索极值并计算纤维方向，步骤如下：
[0028]步骤3.3.1对十二面体进行5次细分，得到10242个球面上的相邻等距点，以球心为 原点得到相应个数的重建向量V，通过3.2)得到的张量系数向量可以求得扩散函数
[0029]
[0030] 步骤3.3.2由前面得到的扩散函数求得在10242个重建向量上的扩散函数值，即纤 维方向分布函数的值，通过搜索纤维方向分布函数值中的极值点来获取纤维的主方向，极 值点的搜索方法如下：
[0031 ]对每一个重建向量vq，q=l，…，10242,在10242个重建向量中搜索出与vq的夹角小 于Θ的所有向量，比较Vq与这些向量所对应的纤维方向分布函数值的大小，若Vq所对应的值 最大，则判断v q为该体素的一个极值方向；依次遍历所有重建向量，最后得到N个极值方向， 这些方向即为纤维的主方向。
[0032] 本发明的有益效果体现在:本发明采用了用高阶笛卡尔张量拟合纤维方向分布函 数平方根的方式，保证了所得纤维方向分布函数的非负性，且角度分辨率较高，实验效果 好。
【具体实施方式】
[0033] 下面对本发明做进一步说明。
[0034] -种非负高阶张量拟牛顿搜索的纤维方向分布估计方法，所述的纤维方向分布估 计方法包括以下步骤：
[0035] (1)数据预处理:读取脑部扩散加权磁共振数据，获取施加梯度方向g时的磁共振 信号S(g)和未施加梯度方向时的磁共振信号So,以及相应的梯度方向数据，选取感兴趣区 域，并计算该区域的扩散衰减信号S(g)/So;
[0036] (2)将感兴趣区域内每个体素中的扩散衰减信号逐个建模为具有扩散形态的椭球 分布模型，建模过程如下：
[0037] 2.1)体素微结构建模:将扩散衰减信号假设为沿重建向量v的单条纤维信号响应 函数R (v，g)与扩散函数D (v)在球面上的卷积：
[0038]
[0039] 其中，j?:卜，@ '「近似看作一个高斯分布函数，g= {gi￡Rlx31 i = l，. . .，n}为梯 度方向，v={vPeRlx3|p = l，. . .，K}为在单位球面上采样的重建向量，R1X3表示维度为1X3 的实数域矩阵空间，η和K分别表示梯度方向和重建向量的个数，μ = ε b是表征扩散效率ε与 扩散敏感系数b共同影响的一个参数，扩散函数D(v)的表达式如下：
[0040]
[0041 ] drs表示单项式v; V, vf Μ的系数，1为高阶张量的阶数，r，s分别表示重建向量v = (vx，Vy，vz)的基方向vx和Vy的指数，λ」表示第j个张量的张量系数，j = l，. . .，m，m表示张量的 个数^(/ +彳+ 2)- ./>) = v:v:v广表示第j个张量单项式，并满足/ = f-I^ + s + 1，: c是由m 个张量系数组成的系数向量，表达式为= 由m个张量单项式构成；
[0042] 2.2)数学模型：
[0043]扩散加权磁共振信号有η个扩散梯度方向gi，i = l，...，n，并且沿重建向量v进行 重建，那么系数向量c通过最小化下面的代价函数J(c)求得：
[0044]
[0045] 其中，Ei = S(gi)/S〇是第i个扩散梯度方向gi上的衰减信号；约>)> 是一个mXm维的矩阵，其值只与扩散梯度方向gl、重建向量v以及参数μ有关，对于每一个扩 散梯度方向gi，都有一个Qi矩阵与之对应；
[0046] (3)计算张量系数向量c，得到扩散函数D(v)，再计算每个采样点处的扩散函数值， 最后将扩散函数值拟合成扩散模型，搜索极值并计算纤维方向。<