本发明涉及基于RL算法的混合响应核函数解决纤维成像中的部分容积效应，利用扩散加权磁共振成像(Diffusion Weighted Magnetic Resonance Imaging，DW-MRI)数据和多响应核函数方法结合RL算法进行纤维方向分布函数的稀疏拟合，从而得到纤维的方向，及在灰质和白质部分容积效应明显的交界处得到更准确的纤维方向。使得到的纤维方向更利于纤维的跟踪。属于医学成像、神经解剖学领域。
背景技术：
核磁共振(MRI)是一种广泛应用在医学成像中的无扩散性方法，作为唯一的活体非入侵式方法，它在帮助人们获得临床神经结构信息和了解大脑皮层区域之间的功能和联系等方面发挥了巨大的作用。脑白质纤维的走向与精神类疾病及脑外科医学疾病存在密切联系，这些信息对于脑的发育、精神分裂症、先天性与获得性脑白质病以及痴呆等的研究提供了新的应用前景。基于扩散加权磁共振成像(DW-MRI)的纤维成像算法能从DW-MRI数据中获得纤维方向信息，为临床医学诊断提供依据，为脑部科学研究提供新的方法。
在各类MRI方法中，扩散张量成像(Diffusion Tensor Imaging，DTI)是较为重要的一种，对于目前已知的多种脑部疾病临床诊断，DTI技术都发挥了不可替代的作用。但传统的DTI方法假设体素内仅含一根纤维，因此无法分辨诸如交叉、瓶颈、分散等复杂的纤维结构，而人脑的神经纤维往往存在交叉、分支或融合的复杂情况，使DTI重构的纤维方向变得不确定。
为了克服DTI的固有局限，高角度分辨率扩散磁共振成像(HARDI)技术应运而生。基于HARDI技术的基础之上，提出了多种纤维重构的方法，例如：Q-ball、扩散谱成像(Diffusion Spectrum Imaging，DSI)、球面反卷积(Spherical deconvolution，SD)等。从目前来看，虽然每种方法都很好的解决了复杂白质部分纤维的成像问题，但是大多数HARDI方法并没有解释非白质(灰质和脑脊液)的部分容积效应对纤维成像的影响。本发明正是利用多响应函数的RL算法来解决非白质部分对脑纤维成像的影响。
技术实现要素：
目前该领域尚未有一种真正意义上解决非白质对脑纤维成像影响的数学模型，为了克服现有方法的不足之处，本发明提出一种基于RL算法的利用多响应核函数来处理非白质的部分容积效应的稀疏成像方法，从而使得白质和灰质交界处的纤维的重构出高分辨率低误差的纤维方向，并进一步区分出脑白质和脑灰质区域。
本发明所采用的技术方案如下：
一种脑纤维稀疏重建的方法，其特征在于：所述方法包括以下步骤：
(1)读取脑部磁共振数据，获取施加梯度方向为g的磁共振扩散信号S(g)，未施加梯度方向的磁共振扩散信号S0及梯度方向数据，对采集的数据进行预处理，选取所需的感兴趣区域，并计算该区域的扩散衰减信号S(g)/S0；
(2)利用Richardson-Lucy迭代算法将感兴趣区域内的每个体素的扩散衰减信号S(g)/S0逐个建模为具有椭球形分布的模型，并增加l1范数正则化对脑纤维进行稀疏重建，建模过程如下：
2.1)体素微结构建模：将扩散衰减信号S(g)/S0假设为沿重建向量v的信号响应核函数H(v,g)与纤维方向分布函数F(v)在球面S2上的卷积：
其中，H(v,g)代表混合响应核函数，它利用脑白质的单条纤维响应核函数和脑灰质中的各向同性响应核函数组成，g＝{gi∈R1×3|i＝1,...,n}为扩散梯度方向向量，v＝{vj∈R1×3|j＝1,...,m}为重建向量，n和m分别为扩散梯度方向向量和重建向量的个数，R是实数集，其数学模型为：
H(v,g)＝faniHani+fisotHisot
其中，,fani，fisot分别是脑白质组织和脑灰质组织的体积分数；分别表示体素中各向异性响应核函数和各向同性响应核函数，各向异性响应核函数Hani是由沿着m个重建方向v的响应核组成，每个响应核都是相同的圆饼状，只是他们的分布方向不同；而各向同性响应核函数也是由沿着m个重建方向v的响应核组成，但其每个响应核的形状都是圆球状；b为扩散敏感系数；表示扩散沿一个主方向进行，在各个方向其扩散程度一致，其中α、β代表纤维扩散程度；
2.2)基于Richardson-Lucy迭代算法的数学模型：
扩散加权磁共振信号有n个扩散梯度方向，并且沿着m个重建向量进行重建，则其数学模型为：
其中，k为迭代次数，F(v)(k)是在当前体素的第k次迭代得到的长度为m×1的列向量，表示沿着m个重建方向均匀分布在球面上的纤维方向分布函数，F(v)(k+1)是当前体素的第k次迭代得到纤维方向分布函数，H即为所述的混合响应核函数H(v,g)，S是在当前体素的包含HARDI信号的长度为n×1的列向量，I0和I1分别是第一类零阶和第一类一阶修正贝塞尔函数，σ2是信号S的方差；
2.3)脑纤维的稀疏重建
用一个完备字典基Φ来表示纤维方向分布函数，即：F(v)＝Φ×c；所得到的系数c恰好是稀疏的，在此基础上，得到了新的Richardson-Lucy算法：
其中c(k)是在当前体素的第k次迭代得到的长度为m×1的系数矩阵，c(k+1)是第k+1次迭代得到的系数矩阵；
2.4)基于Richardson-Lucy迭代算法的l1正则化稀疏重建模型如下：
增加了l1稀疏正则化项，其数学模型为：
其中，L1(k)是第k次迭代的l1正则化项，即是长度为m的列向量，其第i行的元素可以用下式计算得到：
其中，是系数矩阵c(k)的第i行向量在第k次迭代时的梯度方向，和分别表示系数矩阵c(k)的第i行向量对x，y和z方向的偏导，是的二范数，λ是正则化参数；
(3)通过迭代计算得到纤维方向分布函数的系数c，纤维方向分布函数的系数c计算方法包括以下步骤：
3.1)在单位球面上均匀采样m个离散的点，以球心为原点获取这m个重建向量v，计算纤维响应核函数H(v,g)的值，得到n×m的轮换矩阵；
3.2)利用模拟数据模拟仿真，设置迭代初值，令c(0)为各项同性的纤维方向分布函数系数，其振幅设置为1，通过实验选定λ值；
3.3)利用无正则项的Richardson-Lucy算法对感兴趣区的体素进行预处理，得到每个体素的纤维方向分布函数，作为正则化Richardson-Lucy算法的初始纤维方向分布函数值；
3.4)设置迭代终止条件：一是迭代次数；一是迭代误差，令迭代误差为：
所以，迭代次数大于最佳迭代次数或者迭代误差NMSE＜ε作为迭代终止条件；
3.5)通过迭代，得到最优纤维方向分布函数系数cz，其是m列的列向量，再利用完备字典基Φ和最优纤维方向分布函数系数cz重构脑纤维方向分布函数F＝Φ×cz；并利用MATLAB仿真拟合纤维方向分布函数F的分布；
3.6)在数学软件MATLAB中对纤维方向分布函数F进行三维成像，并通过搜索纤维方向分布函数值中的极值点来获取纤维的主方向。
进一步，所述步骤(1)中，所述预处理包括高频滤波、空间降噪和去除线性漂移。
本发明的有益效果为：采用最大似然估计与医学图像处理的理念，因此相比传统方法来看，计算速度快，成像角度分辨率高，且计算鲁棒性高。
附图说明
图1是本发明模拟数据结果图。
图2是本发明实际临床数据效果图。
具体实施方式
下面结合附图对本发明做进一步说明。
参照图1和图2，一种脑纤维稀疏重建的方法，所述方法包括以下步骤：
(1)读取脑部磁共振数据，获取施加梯度方向为g的磁共振扩散信号S(g)，未施加梯度方向的磁共振扩散信号S0及梯度方向数据，对采集的数据进行预处理，选取所需的感兴趣区域，并计算该区域的扩散衰减信号S(g)/S0；
(2)利用Richardson-Lucy迭代算法将感兴趣区域内的每个体素的扩散衰减信号S(g)/S0逐个建模为具有椭球形分布的模型，并增加l1范数正则化对脑纤维进行稀疏重建，建模过程如下：
2.1)体素微结构建模：将扩散衰减信号S(g)/S0假设为沿重建向量v的信号响应核函数H(v,g)与纤维方向分布函数F(v)在球面S2上的卷积：
其中，H(v,g)代表混合响应核函数，它利用脑白质的单条纤维响应核函数和脑灰质中的各向同性响应核函数组成，g＝{gi∈R1×3|i＝1,...,n}为扩散梯度方向向量，v＝{vj∈R1×3|j＝1,...,m}为重建向量，n和m分别为扩散梯度方向向量和重建向量的个数，R是实数集，其数学模型为：
H(v,g)＝faniHani+fisotHisot
其中，fani，fisot分别是脑白质组织和脑灰质组织的体积分数；分别表示体素中各向异性响应核函数和各向同性响应核函数，各向异性响应核函数Hani是由沿着m个重建方向v的响应核组成，每个响应核都是相同的圆饼状，只是他们的分布方向不同；而各向同性响应核函数也是由沿着m个重建方向v的响应核组成，但其每个响应核的形状都是圆球状；b为扩散敏感系数；表示扩散沿一个主方向进行，在各个方向其扩散程度一致，其中α、β代表纤维扩散程度；
2.2)基于Richardson-Lucy迭代算法的数学模型：
扩散加权磁共振信号有n个扩散梯度方向，并且沿着m个重建向量进行重建，则其数学模型为：
其中，k为迭代次数，F(v)(k)是在当前体素的第k次迭代得到的长度为m×1的列向量，表示沿着m个重建方向均匀分布在球面上的纤维方向分布函数，F(v)(k+1)是当前体素的第k次迭代得到纤维方向分布函数，H即为所述的混合响应核函数H(v,g)，S是在当前体素的包含HARDI信号的长度为n×1的列向量，I0和I1分别是第一类零阶和第一类一阶修正贝塞尔函数，σ2是信号S的方差；
2.3)脑纤维的稀疏重建
用一个完备字典基Φ来表示纤维方向分布函数，即：F(v)＝Φ×c；所得到的系数c恰好是稀疏的，在此基础上，得到了新的Richardson-Lucy算法：
其中c(k)是在当前体素的第k次迭代得到的长度为m×1的系数矩阵，c(k+1)是第k+1次迭代得到的系数矩阵；
2.4)基于Richardson-Lucy迭代算法的l1正则化稀疏重建模型如下：
增加了l1稀疏正则化项，其数学模型为：
其中，L1(k)是第k次迭代的l1正则化项，即是长度为m的列向量，其第i行的元素可以用下式计算得到：
其中，是系数矩阵c(k)的第i行向量在第k次迭代时的梯度方向，和分别表示系数矩阵c(k)的第i行向量对x，y和z方向的偏导，是的二范数，λ是正则化参数；
(3)通过迭代计算得到纤维方向分布函数的系数c，纤维方向分布函数的系数c计算方法包括以下步骤：
3.1)在单位球面上均匀采样m个离散的点，以球心为原点获取这m个重建向量v，计算纤维响应核函数H(v,g)的值，得到n×m的轮换矩阵；
3.2)利用模拟数据模拟仿真，设置迭代初值，令c(0)为各项同性的纤维方向分布函数系数，其振幅设置为1，通过实验选定λ值；
3.3)利用无正则项的Richardson-Lucy算法对感兴趣区的体素进行预处理，得到每个体素的纤维方向分布函数，作为正则化Richardson-Lucy算法的初始纤维方向分布函数值；
3.4)设置迭代终止条件：一是迭代次数；一是迭代误差，令迭代误差为：
所以，迭代次数大于最佳迭代次数或者迭代误差NMSE＜ε作为迭代终止条件；
3.5)通过迭代，得到最优纤维方向分布函数系数cz，其是m列的列向量，再利用完备字典基Φ和最优纤维方向分布函数系数cz重构脑纤维方向分布函数F＝Φ×cz；并利用MATLAB仿真拟合纤维方向分布函数F的分布；
3.6)在数学软件MATLAB中对纤维方向分布函数F进行三维成像，并通过搜索纤维方向分布函数值中的极值点来获取纤维的主方向。
图1是本发明模拟数据结果图。其中，模拟数据由下式产生：
其中fi表示第i根纤维所占的比例，f1＝0.5，f2＝0.5，S0＝1，b＝3000s/mm2，扩散张量D的特征值为：λ1＝1.8×10-3mm2/s，λ2＝0.3×10-3mm2/s,λ3＝0.3×10-3mm2/s。81个在半球面内均匀分布的扩散加权磁共振扫描方向，半球面采样点数为321个，图中第一行表示角度，第二行表示纤维方向，第三行表示成像模型，黑线示意两条纤维的方向(通过计算扩散峰值得到)。
图2是本发明实际临床数据效果图。实际数据来自哈佛大学医学院附属医院(Brigham and Women’s Hospital,Brockton VA Hospital,McLean Hospital),利用3-T GE系统和双回波平面成像序列从真实人脑中提取的脑部数据,数据采集参数为:TR＝17000ms,TE＝78ms。体素量为144×144×85,成像域为24cm.平行于AC-PC线的85个轴向切片,每层厚度1.7mm.从51个不同梯度方向扫描数据,扩散敏感系数b＝900s/mm2,8个b＝0的扫描数据。