CN104392019A - 用于脑白质纤维跟踪的高阶扩散张量混合稀疏成像方法 - Google Patents

Abstract

用于脑白质纤维跟踪的高阶扩散张量混合稀疏成像方法，包括以下步骤：读取脑部磁共振数据，获取施加梯度方向g的磁共振信号S(g)，未施加梯度方向的磁共振信号S₀及梯度方向数据，选取所需的感兴趣区域，并计算该区域的扩散衰减信号S(g)/S₀；将感兴趣区域内的每个体素的扩散衰减信号S(g)/S₀逐个建模为具有扩散形态的椭球分布模型；通过计算张量系数λ_j得到扩散函数D(v)，再计算每个采样点的扩散函数值，最后将扩散函数值拟合成扩散模型。本发明涉及压缩感知与稀疏成像的理念，相比传统方法，计算速度快，成像角度分辨率高，且计算样本数可大大减少，实验效果好。

Description

用于脑白质纤维跟踪的高阶扩散张量混合稀疏成像方法

(一)技术领域

本发明涉及图像处理、医学成像、计算方法、数学、三维重建、神经解剖学等领域，尤其是一种用于脑白质纤维跟踪的高阶扩散张量混合稀疏成像方法。

(二)背景技术

脑白质纤维跟踪是获得大脑白质区域纤维走向的一类信息医学技术，它通过追踪局部张量走向来估计纤维的可能路径。就目前而言，脑白质纤维跟踪方法是唯一以种能在活体中无创地获得脑白质纤维走向的方法。纤维跟踪技术首先对原始的DW-MRI数据进行体素建模，获得每个体素内的纤维走向，形成具有解剖学意义的纤维空间微结构，然后再利用纤维跟踪算法对指定区域的纤维方向进行连接。随着磁共振扩散信号的采样精度和对纤维跟踪精度需求的提高，纤维方向估计问题的求解规模越来越大，使得稳定获得高分辨率纤维识别比较困难，这极大地阻碍了该技术在临床医学中的应用。

(三)发明内容

为了克服现有成像方法角度分辨率和计算效率低的缺点，本发明提出一种以高阶张量为导向的高角度分辨率、高效率、低样本数的用于脑白质纤维跟踪的高阶扩散张量混合稀疏成像方法。

本发明所采用的技术方案如下：

用于脑白质纤维跟踪的高阶扩散张量混合稀疏成像方法，所述的混合稀疏成像方法包括以下步骤：

(1)读取脑部磁共振数据，获取施加梯度方向g的磁共振信号S(g)，未施加梯度方向的磁共振信号S₀及梯度方向数据，选取所需的感兴趣区域，并计算该区域的扩散衰减信号S(g)/S₀；

(2)将感兴趣区域内的每个体素的扩散衰减信号S(g)/S₀逐个建模为具有扩散形态的椭球分布模型，建模过程如下：

2.1)体素微结构建模方案：将扩散衰减信号S(g)/S₀假设为沿重建向量v的单条纤维信号响应函数R(v,g)与扩散函数D(v)在球面S²上的卷积：

$S\left(g\right)/S_0=R(v,g)\⊗D\left(v\right)=f_{S^2}R(v,g)D\left(v\right)\mathrm{dv}$

其中，μ表征的是受扩散效率与扩散敏感系数b共同影响的一个参数；表示扩散函数，d_rs表示单项式的系数，λ_j表示第j个张量的张量系数，j＝1,...,m，m表示张量的个数，f_j(v)表示第j个单项式，并满足其中r,s分别表示重建向量v的基方向v₁,v₂的指数；

2.2)构建的数学模型如下：

由于磁共振信号的数据采集并不理想，往往带有一定的噪声，为了尽可能地克服噪声对方向重建的影响，通常利用能量函数的最小化来减小误差；如果扩散加权磁共振信号有n个扩散梯度方向g_i，i＝1,...,n，并且沿重建向量v进行重建，那么λ_j可以通过最小化下面的能量函数E求得

$E=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(S\left(g_i\right)/S_0-\underset{S_2}{\∫}R(v,g_i)D\left(v\right)\mathrm{dv})}^2=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(S\left(g_i\right)/S_0-\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_j\underset{S_2}{\∫}{e}^{-\mathrm{\μ}{\left({g_i}^{T}v\right)}^2}{f}_j\left(v\right)\mathrm{dv})}^2$

为了求解能量函数E得到张量系数λ_j，因此将能量函数方程转换为线性问题

y＝Af

A是n×m的矩阵，其中每个元素张量系数向量f＝[λ₁,λ₂,..λ_m]^T，y是n维列向量，每个元素y_i＝S(g_i)/S₀；

(3)通过计算张量系数λ_j得到扩散函数D(v)，再计算每个采样点的扩散函数值，最后将扩散函数值拟合成扩散模型；张量系数λ_j的计算方法包括以下步骤：

3.1)在单位半球面上均匀采样321个离散的点，以球心为原点获取这321个重建向量v，计算单条纤维响应函数R(v,g)的值，设定一个较低的阶数l_low计算单项式矩阵f(v)；计算R与f(v)的卷积作为步骤2.2)中的A；

3.2)通过稀疏变换过程可以将非稀疏的解映射到稀疏域中，使其稀疏表示，变换过程表示为：

f＝Ψx

其中Ψ为变换矩阵；令Φ＝AΨ，Φ称为传感矩阵；因此我们将上述问题改写为：

y＝Af＝AΨx＝Φx；

3.3)混合加权稀疏算法求解3.2)所提出的问题，其过程包括以下步骤：

步骤3.3.1，在两组线性非负空间-Φ⁰x≥0和x≥0内,通过求解下述约束优化问题获得初始解x⁽⁰⁾来初始化搜索空间：

$\begin{array}{ccccc}{x}^{\left(0\right)}\←\arg \min {\left|\right|{\mathrm{\Φ}}^{0}x-y\left|\right|}_2^2& s.t& -{\mathrm{\Φ}}^{0}x\≥0& \mathrm{and}& x\≥0;\end{array}$

步骤3.3.2，使用低阶混合加权稀疏化方法训练得到正则化矩阵L，即迭代求解下述问题：

$\begin{array}{ccccc}{x}^{(t+1)}\←\arg \min {\left|\right|{\mathrm{\Φ}}^0{x}^{\left(t\right)}-y\left|\right|}_2^2+\mathrm{\ζ}(\mathrm{\α}{\left|\right|{\mathrm{\ω}}^{\left(t\right)}{x}^{\left(t\right)}\left|\right|}_1+(1-\mathrm{\α}){\left|\right|L^{\left(t\right)}{x}^{\left(t\right)}\left|\right|}_2^2)& s.t& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i^{\left(t\right)}=1& \mathrm{and}& x\≥0\end{array}$

其中Φ⁰表示初始的传感矩阵，ω^(t)表示第t次迭代的加权向量，x^(t)表示第t次迭代的解，ζ为约束参数，α用来控制两项罚函数条件||ω^(t)x^(t)||₁和之间的平衡，其中l₁范数||ω^(t)x^(t)||₁用来促进稀疏解的逼近，而l₂范数用来促进解的平滑性；L^(t)表示第t次迭代的正则化矩阵

$L_{m^{\&prime;},n^{\&prime;}}^{\left(t\right)}=\left\{\begin{array}{cc}{D}_{m^{\&prime;},n^{\&prime;}}^{\left(t\right)}& \left|\mathrm{\&mu;}\right|<\mathrm{\&tau;}{\stackrel{\&OverBar;}{D}}^{\left(t\right)}\\ 0& \left|\mathrm{\&mu;}\right|\&GreaterEqual;\mathrm{\&tau;}{\stackrel{\&OverBar;}{D}}^{\left(t\right)}\end{array}\right.$

其中τ为阈值参数，为第t次迭代时的张量D^(t)的平均值，μ＝D^(t)x^(t)，表示D^(t)中满足阈值条件的点，m′,n′表示矩阵中的位置坐标；

步骤3.3.3，更新加权向量ω^(t)：

${\mathrm{\&omega;}}_i^{(t+1)}=\frac{\mathrm{sgn}(x_i^{\left(t\right)}-\mathrm{\&beta;})}{\left|x_i^{\left(t\right)}\right|+\mathrm{\&delta;}}$

其中，表示第t+1次迭代的加权向量ω的第i个元素，sgn是符号函数，β表示符号参数，当时进行符号变换，δ用来防止当为零时分母出现0，当相邻两次迭代的解满足条件时则终止迭代，否则返回步骤3.3.2；

步骤3.3.4，设定一个较高的阶数l_high，按3.1)和3.2)计算高阶下的传感矩阵Φ¹，通过步骤3.3.2训练得到的正则化矩阵L和高阶产生的传感矩阵来求解最终的系数x_finally：

x_finally＝[(Φ¹)^T(Φ¹)+ζ(L)^T(L)]^-1[(Φ¹)^Ty]；

步骤3.3.5，将最终求得的x_finally按步骤3.2)计算张量系数向量f；

3.4)将张量系数用于拟合扩散分布，获取纤维扩散模型，搜索极值计算纤维方向；其过程包括以下步骤：

步骤3.4.1，对十二面体进行四次分形，获取接近于单位球面采样的2562个采样点，得到2562个重建向量v；通过步骤3.3.5所得到的张量系数向量f就可以得到最终的扩散函数D(v)

$D\left(v\right)=\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_j{f}_j\left(v\right);$

步骤3.4.2，通过粒子群局部极值搜索方法获得q个极值点，在每个极值点的邻域附近搜索t个稀疏点；

步骤3.4.3，通过步骤3.4.1和3.4.2得到的扩散函数分别计算这2562个采样点的扩散函数值也称FOD值，而实际上大多数FOD均为0，实际只需计算t个稀疏点的FOD值；使用数学软件MATLAB仿真拟合出FOD值的分布，通过搜索FOD值中的极值点来获取纤维的主方向。

进一步，所述的步骤2.1)中的扩散函数D(v)使用重建向量v的三个基向量v₁,v₂,v₃建立的高阶多项式作为扩散函数，所有的重建向量所构成的扩散函数是一个高阶单项式矩阵。

本发明的有益效果体现在：本发明涉及压缩感知与稀疏成像的理念，相比传统方法，计算速度快，成像角度分辨率高，且计算样本数可大大减少，实验效果好。

(四)附图说明

图1是本发明在模拟磁共振数据下建模的结果图。其中，模拟数据由下式产生：

$S\left(g\right)=S_0\underset{i=1}{\overset{2}{\mathrm{\&Sigma;}}}{f}_i{S}_0{e}^{-b{g}^T\mathrm{Dg}}$

其中f_i表示第i根纤维所占的比例，f₁＝0.5,f₂＝0.5,S₀＝1，扩散敏感系数b＝3000s/mm²，扩散张量D的特征值为：λ₁＝1.8×10^-3mm²/s，λ₂＝0.3×10^-3mm²/s,λ₃＝0.3×10^-3mm²/s，g为梯度方向，实验中l_low＝8,l_high＝18，81个在半球面内均匀分布的扩散加权磁共振扫描方向，半球面采样点数为321个，图中第一行表示角度，第二行表示纤维方向，第三行表示成像模型，图中相交直线示意两条纤维的方向。

图2是本发明采用实际临床数据后的效果图，图中曲线示意了纤维的大致走向。数据来自哈佛大学医学院附属医院(Brigham and Women’s Hospital,Brockton VA Hospital,McLeanHospital)，是利用3-T GE系统和双回波平面成像序列从真实人脑中提取的脑部数据，数据采集参数为：TR＝17000ms，TE＝78ms，体素量为144×144×85，成像域为24cm，平行于AC-PC线的85个轴向切片，每层厚度1.7mm，从51个不同梯度方向扫描数据，扩散敏感系数b＝900s/mm²，8个b＝0的扫描数据。

(五)具体实施步骤

用于脑白质纤维跟踪的高阶扩散张量混合稀疏成像方法，所述的混合稀疏成像方法包括以下步骤：

(1)读取脑部磁共振数据，获取施加梯度方向g的磁共振信号S(g)，未施加梯度方向的磁共振信号S₀及梯度方向数据，选取所需的感兴趣区域，并计算该区域的扩散衰减信号S(g)/S₀；

(2)将感兴趣区域内的每个体素的扩散衰减信号S(g)/S₀逐个建模为具有扩散形态的椭球分布模型，建模过程如下：

2.1)体素微结构建模方案：将扩散衰减信号S(g)/S₀假设为沿重建向量v的单条纤维信号响应函数R(v,g)与扩散函数D(v)在球面S²上的卷积：

$S\left(g\right)/S_0=R(v,g)\&CircleTimes;D\left(v\right)=f_{S^2}R(v,g)D\left(v\right)\mathrm{dv}$

其中，μ表征的是受扩散效率与扩散敏感系数b共同影响的一个参数；表示扩散函数，d_rs表示单项式的系数，λ_j表示第j个张量的张量系数，j＝1,...,m，m表示张量的个数，f_j(v)表示第j个单项式，并满足其中r,s分别表示重建向量v的基方向v₁,v₂的指数；

2.2)构建的数学模型如下：

由于磁共振信号的数据采集并不理想，往往带有一定的噪声，为了尽可能地克服噪声对方向重建的影响，通常利用能量函数的最小化来减小误差；如果扩散加权磁共振信号有n个扩散梯度方向g_i，i＝1,...,n，并且沿重建向量v进行重建，那么λ_j可以通过最小化下面的能量函数E求得

$E=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(S\left(g_i\right)/S_0-\underset{S_2}{\&Integral;}R(v,g_i)D\left(v\right)\mathrm{dv})}^2=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(S\left(g_i\right)/S_0-\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_j\underset{S_2}{\&Integral;}{e}^{-\mathrm{\&mu;}{\left({g_i}^{T}v\right)}^2}{f}_j\left(v\right)\mathrm{dv})}^2$

为了求解能量函数E得到张量系数λ_j，因此将能量函数方程转换为线性问题

y＝Af

A是n×m的矩阵，其中每个元素张量系数向量f＝[λ₁,λ₂,..λ_m]^T，y是n维列向量，每个元素y_i＝S(g_i)/S₀；

(3)通过计算张量系数λ_j得到扩散函数D(v)，再计算每个采样点的扩散函数值，最后将扩散函数值拟合成扩散模型；张量系数λ_j的计算方法包括以下步骤：

3.1)在单位半球面上均匀采样321个离散的点，以球心为原点获取这321个重建向量v，计算单条纤维响应函数R(v,g)的值，设定一个较低的阶数l_low计算单项式矩阵f(v)；计算R与f(v)的卷积作为步骤2.2)中的A；

3.2)通过稀疏变换过程可以将非稀疏的解映射到稀疏域中，使其稀疏表示，变换过程表示为：

f＝Ψx

其中Ψ为变换矩阵；令Φ＝AΨ，Φ称为传感矩阵；因此我们将上述问题改写为：

y＝Af＝AΨx＝Φx；

3.3)混合加权稀疏算法求解3.2)所提出的问题，其过程包括以下步骤：

步骤3.3.1，在两组线性非负空间-Φ⁰x≥0和x≥0内,通过求解下述约束优化问题获得初始解x⁽⁰⁾来初始化搜索空间：

$\begin{array}{ccccc}{x}^{\left(0\right)}\&LeftArrow;\arg \min {\left|\right|{\mathrm{\&Phi;}}^{0}x-y\left|\right|}_2^2& s.t& -{\mathrm{\&Phi;}}^{0}x\&GreaterEqual;0& \mathrm{and}& x\&GreaterEqual;0;\end{array}$

步骤3.3.2，使用低阶混合加权稀疏化方法训练得到正则化矩阵L，即迭代求解下述问题：

$\begin{array}{ccccc}{x}^{(t+1)}\&LeftArrow;\arg \min {\left|\right|{\mathrm{\&Phi;}}^0{x}^{\left(t\right)}-y\left|\right|}_2^2+\mathrm{\&zeta;}(\mathrm{\&alpha;}{\left|\right|{\mathrm{\&omega;}}^{\left(t\right)}{x}^{\left(t\right)}\left|\right|}_1+(1-\mathrm{\&alpha;}){\left|\right|L^{\left(t\right)}{x}^{\left(t\right)}\left|\right|}_2^2)& s.t& \underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i^{\left(t\right)}=1& \mathrm{and}& x\&GreaterEqual;0\end{array}$

其中Φ⁰表示初始的传感矩阵，ω^(t)表示第t次迭代的加权向量，x^(t)表示第t次迭代的解，ζ为约束参数，α用来控制两项罚函数条件||ω^(t)x^(t)||₁和之间的平衡，其中l₁范数||ω^(t)x^(t)||₁用来促进稀疏解的逼近，而l₂范数用来促进解的平滑性；L^(t)表示第t次迭代的正则化矩阵

$L_{m^{\&prime;},n^{\&prime;}}^{\left(t\right)}=\left\{\begin{array}{cc}{D}_{m^{\&prime;},n^{\&prime;}}^{\left(t\right)}& \left|\mathrm{\&mu;}\right|<\mathrm{\&tau;}{\stackrel{\&OverBar;}{D}}^{\left(t\right)}\\ 0& \left|\mathrm{\&mu;}\right|\&GreaterEqual;\mathrm{\&tau;}{\stackrel{\&OverBar;}{D}}^{\left(t\right)}\end{array}\right.$

其中τ为阈值参数，为第t次迭代时的张量D^(t)的平均值，μ＝D^(t)x^(t)，表示D^(t)中满足阈值条件的点，m′,n′表示矩阵中的位置坐标；

步骤3.3.3，更新加权向量ω^(t)：

${\mathrm{\&omega;}}_i^{(t+1)}=\frac{\mathrm{sgn}(x_i^{\left(t\right)}-\mathrm{\&beta;})}{\left|x_i^{\left(t\right)}\right|+\mathrm{\&delta;}}$

其中，表示第t+1次迭代的加权向量ω的第i个元素，sgn是符号函数，β表示符号参数，当时进行符号变换，δ用来防止当为零时分母出现0，当相邻两次迭代的解满足条件时则终止迭代，否则返回步骤3.3.2；

步骤3.3.4，设定一个较高的阶数l_high，按3.1)和3.2)计算高阶下的传感矩阵Φ¹，通过步骤3.3.2训练得到的正则化矩阵L和高阶产生的传感矩阵来求解最终的系数x_finally：

x_finally＝[(Φ¹)^T(Φ¹)+ζ(L)^T(L)]^-1[(Φ¹)^Ty]；

步骤3.3.5，将最终求得的x_finally按步骤3.2)计算张量系数向量f；

3.4)将张量系数用于拟合扩散分布，获取纤维扩散模型，搜索极值计算纤维方向；其过程包括以下步骤：

步骤3.4.1，对十二面体进行四次分形，获取接近于单位球面采样的2562个采样点，得到2562个重建向量v；通过步骤3.3.5所得到的张量系数向量f就可以得到最终的扩散函数D(v)

$D\left(v\right)=\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_j{f}_j\left(v\right);$

步骤3.4.2，通过粒子群局部极值搜索方法获得q个极值点，在每个极值点的邻域附近搜索t个稀疏点；

步骤3.4.3，通过步骤3.4.1和3.4.2得到的扩散函数分别计算这2562个采样点的扩散函数值也称FOD值，而实际上大多数FOD均为0，实际只需计算t个稀疏点的FOD值；使用数学软件MATLAB仿真拟合出FOD值的分布，通过搜索FOD值中的极值点来获取纤维的主方向。

进一步，所述的步骤2.1)中的扩散函数D(v)使用重建向量v的三个基向量v₁,v₂,v₃建立的高阶多项式作为扩散函数，所有的重建向量所构成的扩散函数是一个高阶单项式矩阵。

粒子群搜索方法是进化算法的一种，从随机解出发，通过迭代寻找最优解，具有较高的极值收敛速度。在单位球面采样2562个向量后，对应于扩散函数D(v)可以计算得到2562个函数值，利用粒子群方法在局部搜索空间内搜索尽可能多的极值点，这些极值点所对应的是纤维的主方向。

Claims (2)

1.用于脑白质纤维跟踪的高阶扩散张量混合稀疏成像方法，其特征在于：所述的混合稀疏成像方法包括以下步骤：

(1)读取脑部磁共振数据，获取施加梯度方向g的磁共振信号S(g)，未施加梯度方向的磁共振信号S₀及梯度方向数据，选取所需的感兴趣区域，并计算该区域的扩散衰减信号S(g)/S₀；

(2)将感兴趣区域内的每个体素的扩散衰减信号S(g)/S₀逐个建模为具有扩散形态的椭球分布模型，建模过程如下：

2.1)体素微结构建模方案：将扩散衰减信号S(g)/S₀假设为沿重建向量v的单条纤维信号响应函数R(v,g)与扩散函数D(v)在球面S²上的卷积：

$S\left(g\right)/S_0=R(v,g)\&CircleTimes;D\left(v\right)={\&Integral;}_{S^2}R(v,g)D\left(v\right)\mathrm{dv}$

其中，μ表征的是受扩散效率与扩散敏感系数b共同影响的一个参数；表示扩散函数，d_rs表示单项式的系数，λ_j表示第j个张量的张量系数，j＝1,...,m，m表示张量的个数，f_j(v)表示第j个单项式，并满足其中r,s分别表示重建向量v的基方向v₁,v₂的指数；

2.2)构建的数学模型如下：

由于磁共振信号的数据采集并不理想，往往带有一定的噪声，为了尽可能地克服噪声对方向重建的影响，通常利用能量函数的最小化来减小误差；如果扩散加权磁共振信号有n个扩散梯度方向g_i，i＝1,...,n，并且沿重建向量v进行重建，那么λ_j可以通过最小化下面的能量函数E求得

$E=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(S\left(g_i\right)/S_0-\underset{S_2}{\&Integral;}R(v,g_i)D\left(v\right)\mathrm{dv})}^2=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(S\left(g_i\right)/S_0-\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_j\underset{S_2}{\&Integral;}{e}^{-\mathrm{\&mu;}{\left({g_i}^{T}v\right)}^2}{f}_j\left(v\right)\mathrm{dv})}^2$

为了求解能量函数E得到张量系数λ_j，因此将能量函数方程转换为线性问题

y＝Af

A是n×m的矩阵，其中每个元素张量系数向量f＝[λ₁,λ₂,..λ_m]^T，y是n维列向量，每个元素y_i＝S(g_i)/S₀；

(3)通过计算张量系数λ_j得到扩散函数D(v)，再计算每个采样点的扩散函数值，最后将扩散函数值拟合成扩散模型；张量系数λ_j的计算方法包括以下步骤：

3.1)在单位半球面上均匀采样321个离散的点，以球心为原点获取这321个重建向量v，计算单条纤维响应函数R(v,g)的值，设定一个较低的阶数l_low计算单项式矩阵f(v)；计算R与f(v)的卷积作为步骤2.2)中的A；

3.2)通过稀疏变换过程可以将非稀疏的解映射到稀疏域中，使其稀疏表示，变换过程表示为：

f＝Ψx

其中Ψ为变换矩阵；令Φ＝AΨ，Φ称为传感矩阵；因此我们将上述问题改写为：

y＝Af＝AΨx＝Φx；

3.3)混合加权稀疏算法求解3.2)所提出的问题，其过程包括以下步骤：

步骤3.3.1，在两组线性非负空间-Φ⁰x≥0和x≥0内,通过求解下述约束优化问题获得初始解x⁽⁰⁾来初始化搜索空间：

$\begin{array}{cc}{x}^{\left(0\right)}\&LeftArrow;\arg \min {\left|\right|{\mathrm{\&Phi;}}^{0}x-y\left|\right|}_2^2& s.t-{\mathrm{\&Phi;}}^{0}x\&GreaterEqual;0\mathrm{and\; x}\&GreaterEqual;0;\end{array}$

步骤3.3.2，使用低阶混合加权稀疏化方法训练得到正则化矩阵L，即迭代求解下述问题：

其中Φ⁰表示初始的传感矩阵，ω^(t)表示第t次迭代的加权向量，x^(t)表示第t次迭代的解，为约束参数，α用来控制两项罚函数条件之间的平衡，其中l₁范数||ω^(t)x^(t)||₁用来促进稀疏解的逼近，而l₂范数用来促进解的平滑性；L^(t)表示第t次迭代的正则化矩阵

$L_{m^{\&prime;},n^{\&prime;}}^{\left(t\right)}=\left\{\begin{array}{cc}{D}_{m^{\&prime;},n^{\&prime;}}^{\left(t\right)}& \left|\mathrm{\&mu;}\right|<\mathrm{\&tau;}{\stackrel{\&OverBar;}{D}}^{\left(t\right)}\\ 0& \left|\mathrm{\&mu;}\right|\&GreaterEqual;\mathrm{\&tau;}{\stackrel{\&OverBar;}{D}}^{\left(t\right)}\end{array}\right.$

其中τ为阈值参数，为第t次迭代时的张量D^(t)的平均值，μ＝D^(t)x^(t)，表示D^(t)中满足阈值条件的点，m′,n′表示矩阵中的位置坐标；

步骤3.3.3，更新加权向量ω^(t)：

${\mathrm{\&omega;}}_i^{(t+1)}=\frac{\mathrm{sgn}({x_i}^{\left(t\right)}-\mathrm{\&beta;})}{\left|x_i^{\left(t\right)}\right|+\mathrm{\&delta;}}$

其中，表示第t+1次迭代的加权向量ω的第i个元素，sgn是符号函数，β表示符号参数，当x_i ^(t)＜β时进行符号变换，δ用来防止当x_i ^(t)为零时分母出现0，当相邻两次迭代的解满足条件时则终止迭代，否则返回步骤3.3.2；

步骤3.3.4，设定一个较高的阶数l_high，按3.1)和3.2)计算高阶下的传感矩阵Φ¹，通过步骤3.3.2训练得到的正则化矩阵L和高阶产生的传感矩阵来求解最终的系数x_finally：

步骤3.3.5，将最终求得的x_finally按步骤3.2)计算张量系数向量f；

3.4)将张量系数用于拟合扩散分布，获取纤维扩散模型，搜索极值计算纤维方向；

其过程包括以下步骤：

步骤3.4.1，对十二面体进行四次分形，获取接近于单位球面采样的2562个采样点，得到2562个重建向量v；通过步骤3.3.5所得到的张量系数向量f就可以得到最终的扩散函数D(v)

$D\left(v\right)=\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_j{f}_j\left(v\right);$

步骤3.4.2，通过粒子群局部极值搜索方法获得q个极值点，在每个极值点的邻域附近搜索t个稀疏点；

步骤3.4.3，通过步骤3.4.1和3.4.2得到的扩散函数分别计算这2562个采样点的扩散函数值也称FOD值，而实际上大多数FOD均为0，实际只需计算t个稀疏点的FOD值；使用数学软件MATLAB仿真拟合出FOD值的分布，通过搜索FOD值中的极值点来获取纤维的主方向。

2.如权利要求1所述的方法，其特征在于：所述的步骤2.1)中的扩散函数D(v)使用重建向量v的三个基向量v₁,v₂,v₃建立的高阶多项式作为扩散函数，所有的重建向量所构成的扩散函数是一个高阶单项式矩阵。