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

摘要

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

主权项

用于脑白质纤维跟踪的高阶扩散张量混合稀疏成像方法，其特征在于：所述的混合稀疏成像方法包括以下步骤：(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)\otimes D\left(v\right)={\int }_{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}{\Sigma }}{(S\left(g_i\right)/S_0-\underset{S_2}{\int }R(v,g_i)D\left(v\right)\mathrm{dv})}^2=\underset{i=1}{\overset{n}{\Sigma }}{(S\left(g_i\right)/S_0-\underset{j=1}{\overset{m}{\Sigma }}{\lambda }_j\underset{S_2}{\int }{e}^{-\mu {\left({g_i}^{T}v\right)}^2}{f}_j\left(v\right)\mathrm{dv})}^2$为了求解能量函数E得到张量系数λ_j，因此将能量函数方程转换为线性问题y＝AfA是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|{\Phi }^{0}x-y\left|\right|}_2^2& s.t-{\Phi }^{0}x\ge 0\mathrm{and\; x}\ge 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|\mu \right|<\tau {\bar{D}}^{\left(t\right)}\\ 0& \left|\mu \right|\ge \tau {\bar{D}}^{\left(t\right)}\end{array}\right.$其中τ为阈值参数，为第t次迭代时的张量D^(t)的平均值，μ＝D^(t)x^(t)，表示D^(t)中满足阈值条件的点，m′,n′表示矩阵中的位置坐标；步骤3.3.3，更新加权向量ω^(t)：${\omega }_i^{(t+1)}=\frac{\mathrm{sgn}({x_i}^{\left(t\right)}-\beta )}{\left|x_i^{\left(t\right)}\right|+\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}{\Sigma }}{\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值中的极值点来获取纤维的主方向。