一种脑白质纤维成像的方法
【技术领域】
[0001] 本发明涉及一种脑白质纤维成像的方法，尤其是基于RL算法的混合响应核函 数解决纤维成像中的部分容积效应，利用扩散加权磁共振成像（Diffusion Weighted Magnetic Resonance Imaging，DW-MRI)数据和多响应核函数方法结合RL算法进行纤维方 向分布函数的拟合，从而得到纤维的方向，及在灰质和白质部分容积效应明显的交界处得 到更准确的纤维方向。使得到的纤维方向更利于纤维的跟踪。属于医学成像、神经解剖学 领域。 (二）
【背景技术】
[0002] 核磁共振（MRI)是一种广泛应用在医学成像中的无扩散性方法，作为唯一的活体 非入侵式方法，它在帮助人们获得临床神经结构信息和了解大脑皮层区域之间的功能和联 系等方面发挥了巨大的作用。脑白质纤维的走向与精神类疾病及脑外科医学疾病存在密切 联系，这些信息对于脑的发育、精神分裂症、先天性与获得性脑白质病以及痴呆等的研究提 供了新的应用前景。基于扩散加权磁共振成像（DW-MRI)的纤维成像算法能从DW-MRI数据 中获得纤维方向信息，为临床医学诊断提供依据，为脑部科学研究提供新的方法。
[0003] 在各类MRI方法中，扩散张量成像（Diffusion Tensor Imaging，DTI)是较为重要 的一种，对于目前已知的多种脑部疾病临床诊断，DTI技术都发挥了不可替代的作用。但传 统的DTI方法假设体素内仅含一根纤维，因此无法分辨诸如交叉、瓶颈、分散等复杂的纤维 结构，而人脑的神经纤维往往存在交叉、分支或融合的复杂情况，使DTI重构的纤维方向变 得不确定。
[0004] 为了克服DTI的固有局限，高角度分辨率扩散磁共振成像（HARDI)技术应运而 生。基于HARDI技术的基础之上，提出了多种纤维重构的方法，例如：Q-ball、扩散谱成像 (Diffusion Spectrum Imaging，DSI)、球面反卷积（Spherical deconvolution，SD)等。从 目前来看，虽然每种方法都很好的解决了复杂白质部分纤维的成像问题，但是大多数HARDI 方法并没有解释非白质（灰质和脑脊液）的部分容积效应对纤维成像的影响。目前该领域 尚未有一种真正意义上解决非白质对脑纤维成像影响的数学模型。 (三）
【发明内容】
[0005] 为了克服现有技术中上述的不足之处，本发明提出一种基于RL算法的利用多响 应核函数来处理非白质的部分容积效应的成像方法，从而使得白质和灰质交界处的纤维的 重构出高分辨率低误差的纤维方向，并进一步区分出脑白质和脑灰质区域。
[0006] (1)读取脑部磁共振数据，获取施加梯度方向为g的磁共振扩散信号S (g)，未施 加梯度方向的磁共振扩散信号S。及梯度方向数据，对采集的数据进行预处理，包括：高频 滤波，空间降噪，去除线性漂移等；选取所需的感兴趣区域，并计算该区域的扩散衰减信号 S(g)/S0;
[0007] (2)利用Richardson-Lucy迭代算法将感兴趣区域内的每个体素的扩散衰减信号 S (g) /S。逐个建模为具有椭球形分布的模型，并增加全变差正则化去除噪声影响，建模过程 如下：
[0008] 2. 1)体素微结构建模方案：将扩散衰减信号S (g) /S。假设为沿重建向量V的信号 响应核函数R(V，g)与扩散方向分布函数f (V)在球面S2上的卷积：
[0010] R(v，g)代表混合响应核函数，它利用脑白质的单纤维响应核函数和脑灰质中 的各向同性响应核函数组成，可以有效的解决脑白质纤维成像中的部分容积效应，g = {gie R 1X3| i = 1，...，η}为扩散梯度方向，V = {vf R 1X3| j = 1，…，m}为重建向量；其 数学模型为：
分别表不体素中各向异性响应核函数和各向同性响 应核函数，b是扩散敏感系数；各向异性响应核函数Rani是由沿着m个重建方向V的响应核 组成，每个响应核都是相同的圆饼状，只是他们的分布方向不同；而各向同性响应核函数只 由一个单独的响应核组成，其形状是圆球状；
表示扩散沿一个主方向进行，
在各个方向其扩散程度一致，其中α、β代表纤维扩散程度，其数值可在相 关文献中查阅得到；
[0013] 2. 2) 一种基于Richardson-Lucy迭代算法的数学模型如下：
[0014] 由于磁共振信号的数据采集并不理想，往往带有一定的噪声，为了尽可能地克服 噪声对纤维方向重建的影响，采用经典的图像恢复方法中的Richardson-Lucy迭代算法， 有效的保证重构的纤维方向分布函数的非负性，减少重构伪峰的数量；该数学模型是通过 最大似然估计的方法计算获得；如果扩散加权磁共振信号有η个扩散梯度方向，并且沿着m 个重建向量进行重建，则其数学模型为：
[0016] 其中k为迭代次数，i表示向量的第i行，f(k)是第k次迭代得到的长度为m的列 向量，表不沿着m个重建方向均勾分布在球面上的扩散方向分布函数，S是包含HARDI信号 的长度为η的列向量；
[0017] 2. 3)基于Richardson-Lucy迭代算法的TV正则化模型如下：
[0018] Richardson-Lucy算法本身具有局限性，虽然该算法可以在一定程度上抑制噪声 的对成像的影响；但随着迭代次数的增加，由Richardson-Lucy算法计算得到的最优解会 被噪声所控制；为了克服迭代次数增加带来的不良影响，在此算法的基础之上增加了全变 差正则化项来抑制噪声，其数学模型为：
[0022] TV1(k)是第k次迭代的全变差正则化项；▽/?是第k次迭代的图像扩散方向分布函 数f(k)的梯度，div表示散度，λ是正则化参数；
[0023] (3)通过迭代计算得到扩散方向分布函数f (V)，扩散方向分布函数f (V)的计算方 法包括以下步骤：
[0024] 3. 1)在单位球面上均匀采样m个离散的点，以球心为原点获取这m个重建向量V， 计算纤维响应核函数R(V，g)的值，得到ηX (m+Ι)的轮换矩阵；
[0025] 3. 2)利用模拟数据模拟仿真，设置迭代初值；令fw为各项同性的扩散方向分布 函数，其振幅设置为1 ;因为fW初始值是正的，所以算法的非负性限定自然的满足；可以通 过实验选定λ值，来平衡数据项和正则项对迭代算法的影响；
[0026] 3. 3)利用无正则项的RL算法对感兴趣区的体素进行预处理；得到每个体素的扩 散方向分布函数f，作为正则化RL算法的初始扩散方向分布函数值；
[0027] 3. 4)设置迭代终止条件：一是迭代次数；一是迭代误差，令迭代误差为：
[0029] 所以，迭代次数大于最佳迭代次数或者迭代误差d〈 ε (-般ε取10 3)作为迭代 终止条件；（对于模拟数据，已知最优的扩散方向分布函数解f，所以上式迭代误差中f(k) = f ;)
[0030] 3.5)通过迭代，得到了最优的扩散方向分布函数f，其是m+1列的列向量，其中最 后一列的值4+1)即为每个体素中各向同性扩散所占的相对比例大小；而前m列即构成体素 中各向异性纤维方向分布函数f〇d，即纤维方向分布函数fod是扩散方向分布函数f的前 m-Ι项；通过实验可以知道：当4+1)> Θ时，体素中的扩散呈各向同性；不对其进行纤维方向 的重构，而对于f Θ的体素，利用MATLAB仿真拟合纤维方向分布函数fod的方向；
[0031] 3. 6)在数学软件MATLAB中对纤维方向分布函数fod进行三维成像，并通过搜索纤 维方向分布函数值中的极值点来获取纤维的主方向。
[0032] 本发明正是利用多响应函数的RL算法来解决非白质部分对脑纤维成像的影响。 本发明涉及最大似然估计与医学图像处理的理念；
[0033] 本发明的优点是：计算速度快，成像角度分辨率高，且计算鲁棒性高。据笔者经验 看来，本发明所得效果是目前该领域最好效果。 (四）【附图说明】
[0034] 图1是本发明模拟数据结果图。其中，模拟数据由下式产生：
[0036] 其中f表示第i跟纤维所占的比例，
3000s/mm2,扩散张量 D 的特征值为：A1= 1.8X10 3mm2/s，λ2= 0.3X10 3mm2/s, λ3 = 0. 3Χ 10 3mm2/s。81个在半球面内均勾分布的扩散加权磁共振扫描方向，半球面采样点数为 321个，图中第一行表示角度，第二行表示纤维方向，第三行表示成像模型，黑线示意两条纤 维的方向（通过计算扩散峰值得到）。
[0037] 图2是本发明实际临床数据效果图。实际数据来自哈佛大学医学院附属医院 (Brigham and Women' s Hospital, Brockton VA Hospital, McLean Hospital),利用 3-T GE系统和双回波平面成像序列从真实人脑中提取的脑部数据，数据采集参数为：TR = 17000ms, TE = 78ms。体素量为144 X 144X85,成像域为24cm.平行于AC-PC线的85个轴 向切片，每层厚度I. 7mm.从51个不同梯度方向扫描数据，扩散敏感系数b = 900s/mm2, 8 个b=0的扫描数据。
[0038] (五）具体实施步骤
[0039] 下面结合附图进一步说明本发明。
[0040] 本发明所述的一种脑白质纤维成像的方法，包括如下步骤：
[0041] (1)读取脑部磁共振数据，获取施加梯度方向为g的磁共振扩散信号S (g)，未施 加梯度方向的磁共振扩散信号S。及梯度方向数据，对采集的数据进行预处理，包括：高频 滤波，空间降噪，去除线性漂移等；选取所需的感兴趣区域，并计算该区域的扩散衰减信号 S(g)/S0;
[0042] (2)利用Richardson-Lucy迭代算法将感兴趣区域内的每个体素的扩散衰减信号 S (g) /S。逐个建模为具有椭球形分布的模型，并增加全变差正则化去除噪声影响，建模过程 如下：
[0043] 2. 1)体素微结构建模方案：将扩散衰减信号S (g) /S。假设为沿重建向量V的信号 响应核函数R(V，g)与扩散方向分布函数f (V)在球面S2上的卷积：
[0045] R(v，g)代表混合响应核函数，它利用脑白质的单纤维响应核函数和脑灰质中 的各向同性响应核函数组成，可以有效的解决脑白质纤维成像中的部分容积效应，g = {gie R 1X3| i = 1，...，η