一种脑白质纤维成像的方法

摘要

利用混合响应核函数解决脑白质纤维成像的部分容积效应的方法，包括以下步骤：读取脑部磁共振数据，获取施加梯度方向g的磁共振信号S(g)，未施加梯度方向的磁共振信号S₀及梯度方向数据，选取所需的感兴趣区域，并计算该区域的扩散衰减信号S(g)/S₀；利用Richardson‐Lucy迭代算法将感兴趣区域内的每个体素的扩散衰减信号S(g)/S₀逐个建模为具有椭球形分布的模型，并增加全变差正则化去除噪声影响，通过Richardson‐Lucy迭代算法得到扩散方向分布函数f(v)，从f(v)中分离出各向同性扩散值f_(m+1)和各向异性扩散的纤维分布函数fod。本发明涉及概率论与全变差正则化的理念，相比传统方法，计算速度快，成像角度分辨率高，可以区分出脑灰质和脑白质，实验效果好。

主权项

一种脑白质纤维成像的方法，其成像方法包括以下步骤：(1)读取脑部磁共振数据，获取施加梯度方向为g的磁共振扩散信号S(g)，未施加梯度方向的磁共振扩散信号S₀及梯度方向数据，对采集的数据进行预处理，包括：高频滤波，空间降噪，去除线性漂移等；选取所需的感兴趣区域，并计算该区域的扩散衰减信号S(g)/S₀；(2)利用Richardson‑Lucy迭代算法将感兴趣区域内的每个体素的扩散衰减信号S(g)/S₀逐个建模为具有椭球形分布的模型，并增加全变差正则化去除噪声影响，建模过程如下：2.1)体素微结构建模方案：将扩散衰减信号S(g)/S₀假设为沿重建向量v的信号响应核函数R(v,g)与扩散方向分布函数f(v)在球面S²上的卷积：$S\left(g\right)/S_0=R(v,g)\otimes f\left(v\right)={\int }_{S^2}R(v,g)f\left(v\right)dv$R(v,g)代表混合响应核函数，它利用脑白质的单纤维响应核函数和脑灰质中的各向同性响应核函数组成，可以有效的解决脑白质纤维成像中的部分容积效应，g＝{g_i∈R^1×3|i＝1,...,n}为扩散梯度方向，v＝{v_j∈R^1×3|j＝1,...,m}为重建向量；其数学模型为：R(v,g)＝[R_ani R_isot]其中，分别表示体素中各向异性响应核函数和各向同性响应核函数，b是扩散敏感系数；各向异性响应核函数R_ani是由沿着m个重建方向v的响应核组成，每个响应核都是相同的圆饼状，只是他们的分布方向不同；而各向同性响应核函数只由一个单独的响应核组成，其形状是圆球状；$D_{ani}=\left[\begin{array}{ccc}\alpha & 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$表示扩散沿一个主方向进行，$D_{iso}=\left[\begin{array}{ccc}\beta & 0& 0\\ 0& \beta & 0\\ 0& 0& \beta \end{array}\right]$在各个方向其扩散程度一致，其中α、β代表纤维扩散程度，其数值可在相关文献中查阅得到；2.2)一种基于Richardson‑Lucy迭代算法的数学模型如下：由于磁共振信号的数据采集并不理想，往往带有一定的噪声，为了尽可能地克服噪声对纤维方向重建的影响，采用经典的图像恢复方法中的Richardson‑Lucy迭代算法，有效的保证重构的纤维方向分布函数的非负性，减少重构伪峰的数量；该数学模型是通过最大似然估计的方法计算获得；如果扩散加权磁共振信号有n个扩散梯度方向，并且沿着m个重建向量进行重建，则其数学模型为：$f{{\left(v\right)}_i}^{\left(k+1\right)}=f{{\left(v\right)}_i}^{\left(k\right)}\frac{{\left[R{(v,g)}^{T}S\right]}_i}{{\left[R{(v,g)}^{T}R(v,g)f{\left(v\right)}^{\left(k\right)}\right]}_i}$其中k为迭代次数，i表示向量的第i行，f^(k)是第k次迭代得到的长度为m的列向量，表示沿着m个重建方向均匀分布在球面上的扩散方向分布函数，S是包含HARDI信号的长度为n的列向量；2.3)基于Richardson‑Lucy迭代算法的TV正则化模型如下：Richardson‑Lucy算法本身具有局限性，虽然该算法可以在一定程度上抑制噪声的对成像的影响；但随着迭代次数的增加，由Richardson‑Lucy算法计算得到的最优解会被噪声所控制；为了克服迭代次数增加带来的不良影响，在此算法的基础之上增加了全变差正则化项来抑制噪声，其数学模型为：$f{{\left(v\right)}_i}^{\left(k+1\right)}=f{{\left(v\right)}_i}^{\left(k\right)}\frac{{\left[R{(v,g)}^{T}S\right]}_i}{{\left[R{(v,g)}^{T}R(v,g)f{\left(v\right)}^{\left(k\right)}\right]}_i}\times {{\mathrm{TV}}_i}^{\left(k\right)}$其中，${{\mathrm{TV}}_i}^{\left(k\right)}=\frac{1}{1-\lambda div\left(\frac{▿{\left[f^{\left(k\right)}\right]}_i}{\left[\right|▿f^{\left(k\right)}{|}_i]}\right)}$TV_i^(k)是第k次迭代的全变差正则化项；是第k次迭代的图像扩散方向分布函数f^(k)的梯度，div表示散度，λ是正则化参数；(3)通过迭代计算得到扩散方向分布函数f(v)，扩散方向分布函数f(v)的计算方法包括以下步骤：3.1)在单位球面上均匀采样m个离散的点，以球心为原点获取这m个重建向量v，计算纤维响应核函数R(v,g)的值，得到n×(m+1)的轮换矩阵；3.2)利用模拟数据模拟仿真，设置迭代初值；令f⁽⁰⁾为各项同性的扩散方向分布函数，其振幅设置为1；因为f⁽⁰⁾初始值是正的，所以算法的非负性限定自然的满足；可以通过实验选定λ值，来平衡数据项和正则项对迭代算法的影响；3.3)利用无正则项的RL算法对感兴趣区的体素进行预处理；得到每个体素的扩散方向分布函数f，作为正则化RL算法的初始扩散方向分布函数值；3.4)设置迭代终止条件：一是迭代次数；一是迭代误差，令迭代误差为：$d=\frac{\left|\right|f^{(k+1)}-f^{\left(k\right)}\left|\right|}{\left|\right|f^{\left(k\right)}\left|\right|}$所以，迭代次数大于最佳迭代次数或者迭代误差d<ε(一般ε取10^‑3)作为迭代终止条件；(对于模拟数据，已知最优的扩散方向分布函数解f，所以上式迭代误差中f^(k)＝f；)3.5)通过迭代，得到了最优的扩散方向分布函数f，其是m+1列的列向量，其中最后一列的值f_(m+1)即为每个体素中各向同性扩散所占的相对比例大小；而前m列即构成体素中各向异性纤维方向分布函数fod，即纤维方向分布函数fod是扩散方向分布函数f的前m‑1项；通过实验可以知道：当f_(m+1)>θ时，体素中的扩散呈各向同性；不对其进行纤维方向的重构，而对于f_(m+1)<θ的体素，利用MATLAB仿真拟合纤维方向分布函数fod的方向；3.6)在数学软件MATLAB中对纤维方向分布函数fod进行三维成像，并通过搜索纤维方向分布函数值中的极值点来获取纤维的主方向。