一种基于空间结构一致性的脑纤维微结构重构方法

摘要

一种基于空间结构一致性的脑纤维微结构重构方法，包括如下步骤：步骤一，采用基于球面字典反褶积的双瓣基函数，步骤二，建立空间结构一致性模型，基于可分空间域使用数学软MATLAB仿真拟合出FOD值的分布，通过搜索FOD值中的极值点来获取纤维的主方向，实现脑纤维微结构重构。本发明提供一种提升计算效率、精度较高的基于空间结构一致性的脑纤维微结构重构方法。

主权项

一种基于空间结构一致性的脑纤维微结构重构方法，其特征在于：该重构方法包括如下步骤：步骤一，采用基于球面字典反褶积的双瓣基函数：在旋转向量和中心向量上的纤维概率分布函数f(v|u)称为fODF，其中m_v表示旋转向量组的空间维度、m_u表示的是中心向量组的空间维度，通过球面反褶积方法将纤维解剖结构描述为核函数的卷积，在扩散梯度g∈S²上的测量信号s(g|u)：$s\left(g\right|u)={\int }_{S^2}r(g,v)f\left(v\right|u)d\mu \left(v\right)$其中r(g,v)表示核函数，μ(v)是在S²上的Haar测度；方便起见，s(g|u)信号可由均匀分布在单位球面上的离散采样μ来表示：$s\left(g\right|u)=\underset{i=1}{\overset{m_u}{\Sigma }}r(g,v)f\left(v\right|u_i)$定义纤维响应函数：$r(g,v)=e^{-l{\left({\mathrm{gv}}^T\right)}^2}$其中是表征扩散张量D的对角矩阵，λ_max是该对角矩阵中主要特征向量，l表示扩散敏感系数以及各向异性相互作用对信号衰减的影响程度；以此扩展出在多壳采样方案上的纤维响应函数的估计：$r(g,v)=\{r(g_{b_J},v)=e^{-l_{b_J}{\left(g_{b_J}{v}^T\right)}^2}|J=1,...,q\}$其中J＝1,...,q指代第J个球壳，对应分布在b_J壳上的梯度向量集g，l_bJ表示第J个球壳下的第bJ个扩散敏感系数所指代的l；从而推得一种新的纤维方向分布函数形式：$f\left(v\right|u)=\underset{i=1}{\overset{m_u}{\Sigma }}{w}_{i}d(v,u_i),u_i\in u$d(v,u_i)表示向量组v和第u_i个方向下的一组过完备的方向分布基函数，w_i是基函数的相对权重，是w_i中非零元素的数目；表观扩散系数可以被近似估计：g^TDg≈λ_maxg^T(vv^T)g＝λ_maxcos²θ(g,v)θ(g,v)表示夹角余弦函数，D是表示由所有的d(v,u_i)构成的一组基函数字典；为了描述球面坐标下的字典基分布，我们将原有直角坐标系下的d(v,u_i)描述为球面坐标下的是离散集上的最小角余弦，κ₁将纤维方向分布函数标准化至单位球，κ₂是用于调整峰值的参数，τ表示偶次幂；纤维方向分布函数(fODF)中基函数系数w的估计：$\underset{w}{\arg min}\left|\right|\Phi w-s|{|}_2^2$Φ是观测矩阵，s是信号向量；为避免伪峰与复杂的算法以及更高阶的大规模病态逆问题的计算，直接通过非负最小二乘方法求得基函数系数w的估计：$\underset{w\ge 0}{w^{*}\leftarrow \arg min}\left|\right|\Phi w-s|{|}_2^2$w^*表示w的最优解；步骤二，建立空间结构一致性模型：利用贝叶斯公式，得到：P(x|s)∝P(s|x)P(x)后验概率密度P(x|s)正比于数据似然函数P(s|x)和先验概率密度函数P(x)的乘积；之后改写为：$P\left(x\right|s)\propto 1/(U^{In}+\underset{i}{\Sigma }{\beta }_i{U}_i^{Ex})$U^In是内部能量，U^Ex是外部能量，β是先验分布的超参数；后验概率的最大化转化为总能量函数的最小化：$\arg \max P\left(x\right|s)\stackrel{\Delta }{=}\mathrm{argmin}\{U^{In}+\underset{i}{\Sigma }{\beta }_i{U}_i^{Ex}\}$P(x|s)是后验概率密度，U^In是内部能量，U^Ex是外部能量，β是先验分布的超参数；其中内部能量：$U^{In}\propto \left|\right|S-S^{\prime }|{|}_2^2=||S-\Theta W|{|}_2^2$S是测量信号集，S′表示的是待估计的信号，W是w系数组，Θ是基于观测矩阵的块对角矩阵；外部能量：$U^{Ex}=\underset{c\in \Omega ,\tilde{c}\in H_c⋐\Omega }{\Sigma }\left|\right|M(w_c-{\bar{w}}_{\tilde{c}})|{|}_2$U^Ex表示外部能量，表示纤维方向分布函数的算数平均值，M通过降采样方向v_t去训练一个字典基来获得，w_c表示第c∈Ω个体素的系数；W的估计：$\underset{w_c\ge 0}{\arg min}\left\{\right||s_c-{\mathrm{\Phi w}}_c|{|}_2^2+{\beta }_c\left|\right|M(w_c-{\bar{w}}_{\tilde{c}})\left|{|}_2\right\},c\in \Omega ;\tilde{c}\in H_c⋐\Omega$s_c是测量信号在体素c中的系数，w_c是字典在体素c中的系数；为了得到体素交叉纤维的结构，定义全局代价函数：$\underset{W\ge 0}{\arg \min }\left\{\left|\right|S-\Theta W|{|}_2^2+\underset{c\in \Omega }{\Sigma }\left(\underset{\tilde{c}\in H_c⋐\Omega }{\Sigma }{\sigma }_{c,\tilde{c}}\left|\right|w_c-w_{\tilde{c}}|{|}_2+{\beta }_1|\left|Q\left(w_c-{\tilde{w}}_{\tilde{c}}\right)\right||+J\left\{w_c|c\in \Omega \right\}\right)\right\}$通过计算w_c和w_c周围邻域的系数余弦距离得到，β₁是人为定义一个参数，Q是由在基函数上训练字典得到；空间结构一致性模型的局部线性近似估计：$\begin{array}{c}\left(w_c^{\left(t+1\right)},{\xi }_c^{\left(t+1\right)}\right)=\underset{{\xi }_c^{\left(t+1\right)},w_c^{\left(t+1\right)}\ge 0}{\arg \min }\left\{\right||s_c-{\mathrm{\Phi w}}_c^{\left(t\right)}|{|}_2^2+\underset{\tilde{c}\in H_c⋐\Omega }{\Sigma }{\sigma }_{c,\tilde{c}}\left|\right|w_c^{\left(t\right)}-w_{\tilde{c}}|{|}_2\\ +{\beta }_1\left|\right|Q\left(w_c^{\left(t\right)}-{\tilde{w}}_{\tilde{c}}\right)|{|}_2+{\beta }_2|\left|{\xi }_c^{\left(t\right)}\right|{|}_1+{\delta }_{\xi }\left|\right|w_c^{\left(t\right)}-{\xi }_c^{\left(t\right)}+p_{c,\xi }^{\left(t\right)}\left|{|}_2\right\}\\ p_{c,\xi }^{\left(t+1\right)}=p_{c,\xi }^{\left(t\right)}+w_c^{\left(t+1\right)}-{\xi }_c^{\left(t+1\right)}\end{array}$t是迭代指数，δ_ξ是预定义辅助常数，是第t次迭代的扩展拉格朗日乘子向量，上式分拆为两部分：$\begin{array}{c}{w}_c^{\left(t+1\right)}=\underset{w_c^{\left(t+1\right)}\ge 0}{\arg \min }\left\{\right||s_c-{\mathrm{\Phi w}}_c^{\left(t\right)}|{|}_2^2+\underset{\tilde{c}\in H_c⋐\Omega }{\Sigma }{\sigma }_{c,\tilde{c}}\left|\right|w_c^{\left(t\right)}-w_{\tilde{c}}|{|}_2\\ +{\beta }_1\left|\right|Q\left(w_c^{\left(t\right)}-{\tilde{w}}_{\tilde{c}}\right)|{|}_2+{\delta }_{\xi }||w_c^{\left(t\right)}-{\xi }_c^{\left(t\right)}+p_{c,\xi }^{\left(t\right)}|{|}_2\}\end{array}$${\xi }_c^{(t+1)}=\underset{{\xi }_c^{\left(t+1\right)}\ge 0}{\mathrm{argmin}}\left\{{\delta }_{\xi }\right||w_c^{\left(t\right)}-{\xi }_c^{\left(t\right)}+p_{c,\xi }^{\left(t\right)}|{|}_2+{\beta }_2\left|\right|{\xi }_c^{\left(t\right)}\left|{|}_1\right\}$优化为一个可分空间域，用增强拉格朗日方法解决，得到：$\begin{array}{c}{w}_c^{\left(t+1\right)}={\left({\Phi }^T\Phi +\underset{\tilde{c}\in H_c⋐\Omega }{\Sigma }{\sigma }_{c,\tilde{c}}I+{\beta }_c{Q}^{T}Q+{\delta }_{\xi }I\right)}^{-1}·\\ {\left({\Phi }^T{s}_c+\underset{\tilde{c}\in H_c⋐\Omega }{\Sigma }{\sigma }_{c,\tilde{c}}{w}_{\tilde{c}}+{\beta }_1{Q}^{T}Q{\tilde{w}}_{\tilde{c}}+{\delta }_{\xi }\left({\xi }_c^{\left(t\right)}-p_{c,\xi }^{\left(t\right)}\right)I\right)}_{w_c^{\left(t+1\right)}\ge 0}\end{array}$其中I表示单位矩阵；基于可分空间域使用数学软MATLAB仿真拟合出FOD值的分布，通过搜索FOD值中的极值点来获取纤维的主方向。