[发明专利]一种非负高阶张量拟牛顿搜索的纤维方向分布估计方法

权利要求书

1.一种非负高阶张量拟牛顿搜索的纤维方向分布估计方法，其特征在于：所述的纤维方向分布估计方法包括以下步骤：

(1)数据预处理：读取脑部扩散加权磁共振数据，获取施加梯度方向g时的磁共振信号S(g)和未施加梯度方向时的磁共振信号S0，以及相应的梯度方向数据，选取感兴趣区域，并计算该区域的扩散衰减信号S(g)/S0；

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

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

S ( g ) / S 0 = R ( v , g ) ⊗ D ( v ) = ∫ R ( v , g ) D ( v ) d v ]]>

其中，近似看作一个高斯分布函数，g＝{gi∈R1×3|i＝1,...,n}为梯度方向，v＝{vp∈R1×3|p＝1,...,K}为在单位球面上采样的重建向量，R1×3表示维度为1×3的实数域矩阵空间，n和K分别表示梯度方向和重建向量的个数，μ＝εb是表征扩散效率ε与扩散敏感系数b共同影响的一个参数，扩散函数D(v)的表达式如下：

D ( v ) = ( Σ r = 0 l Σ s = 0 l - r d r s v x r v y s v z l - r - s ) 2 = ( Σ j = 1 m λ j f j ( v ) ) 2 = ( F ( v ) c ) 2 = c T ( F T ( v ) F ( v ) ) c ]]>

drs表示单项式的系数，l为高阶张量的阶数，r,s分别表示重建向量v＝(vx,vy,vz)的基方向vx和vy的指数，λj表示第j个张量的张量系数，j＝1,...,m，m表示张量的个数表示第j个张量单项式，并满足c是由m个张量系数组成的系数向量，表达式为c＝[λ1,λ2,..λm]T，F(v)＝[f1(v),f2(v),...fm(v)]由m个张量单项式构成；

2.2)数学模型：

扩散加权磁共振信号有n个扩散梯度方向gi，i＝1,...,n，并且沿重建向量v进行重建，那么系数向量c通过最小化下面的代价函数J(c)求得：

J ( c ) = Σ i = 1 n ( ∫ R ( v , g i ) D ( v ) d v - S ( g i ) / S 0 ) 2 = Σ i = 1 n ( ∫ e - μ ( g i T v ) 2 ( c T F T ( v ) F ( v ) c ) d v - S ( g i ) / S 0 ) 2 = Σ i = 1 n ( c T ( ∫ e - μ ( g i T v ) 2 ( F T ( v ) F ( v ) ) d v ) c - S ( g i ) / S 0 ) 2 = Σ i = 1 n ( cQ i c T - E i ) 2 ]]>

其中，Ei＝S(gi)/S0是第i个扩散梯度方向gi上的衰减信号；是一个m×m维的矩阵，其值只与扩散梯度方向gi、重建向量v以及参数μ有关，对于每一个扩散梯度方向gi，都有一个Qi矩阵与之对应；

(3)计算张量系数向量c，得到扩散函数D(v)，再计算每个采样点处的扩散函数值，最后将扩散函数值拟合成扩散模型，搜索极值并计算纤维方向；

所述步骤(3)中，所述张量系数向量c的计算包括以下步骤：

3.1)在单位半球面上均匀采样321个离散的点，以球心为原点获取这321个重建向量v，计算单条纤维响应函数R(v,g)的值，设定高阶张量模型的阶数l，计算单项式矩阵F(v)，进而计算出步骤2.2)中的矩阵Qi；

3.2)使用BFGS拟牛顿搜索算法迭代求解2.2)中的最小化问题，步骤如下：

步骤3.2.1已知代价函数J(c)，选取一个初始点c1作为第一次迭代的搜索起始点，迭代次数计为k＝1，并设置最大迭代次数k_max，计算J(c)的梯度向量表达式：

∂ J ( c k ) ∂ c k = 4 Σ i = 1 N ( c k T Q i c k - E i ) Q i c k ]]>

步骤3.2.2对于第k次迭代，令表示第k次迭代的梯度向量，计算拟牛顿方向dk＝-Hktk作为第k次迭代的搜索方向，以ck为起点，沿方向dk进行一维搜索，求得本次搜索的可接受步长αk，矩阵Hk表示的是第k次迭代时ck点处Hesse矩阵的逆的近似矩阵，其更新方法如下：

其中，I是单位阵，δk＝ck+1-ck是相邻两次迭代的解的差，βk＝tk+1-tk是相邻两次迭代的梯度向量之差；

步骤3.2.3更新张量系数向量ck+1＝ck+αkdk，计算对应的代价函数值J(ck+1)，当达到最大迭代次数或相邻两次迭代的代价函数值满足时，则终止迭代，否则，返回步骤3.2.2，进入下一次迭代；参数σ为一个较小的常数，用于判断算法是否收敛到一个局部极小值。

2.如权利要求1所述的一种非负高阶张量拟牛顿搜索的纤维方向分布估计方法，其特征在于：所述步骤(3)还包括3.3)将得到的张量系数用于拟合扩散函数，获取纤维方向分布函数模型，搜索极值并计算纤维方向，步骤如下：

步骤3.3.1对十二面体进行5次细分，得到10242个球面上的相邻等距点，以球心为原点得到相应个数的重建向量V，通过3.2)得到的张量系数向量可以求得扩散函数

D ( V ) = ( Σ j = 1 m λ j f j ( V ) ) 2 ]]>

步骤3.3.2由前面得到的扩散函数求得在10242个重建向量上的扩散函数值，即纤维方向分布函数的值，通过搜索纤维方向分布函数值中的极值点来获取纤维的主方向，极值点的搜索方法如下：

对每一个重建向量Vq，q＝1,…,10242，在10242个重建向量中搜索出与Vq的夹角小于θ的所有向量，比较Vq与这些向量所对应的纤维方向分布函数值的大小，若Vq所对应的值最大，则判断Vq为体素的一个极值方向；依次遍历所有重建向量，最后得到N个极值方向，这些方向即为纤维的主方向。