基于非局部均值平滑的体素内纤维结构估计

楚春雨

渤海大学工学院，辽宁锦州121013

前言

磁共振扩散成像是目前唯一能够在活体上测量组织内水分子扩散运动与成像的无创方法，它通过测量和量化组织中水分子的扩散信息来探测组织的微观结构。水分子沿不同方向的扩散信息包含在一组不同扩散梯度方向的扩散加权图像（Diffusion Weighted Image,DWI）中，通过对扩散函数进行建模可以解析出每个体素内的纤维结构信息（主要是纤维的走向）。根据每个体素内的纤维走向，利用纤维束追踪技术可重建出组织纤维束的整体三维结构，并可从中提取有效的临床统计特征，从而用于医疗诊断及相关研究等。

扩散张量成像（Diffusion Tensor Imaging,DTI）［1］是最早提出的一种扩散模型，目前已被广泛用于临床和医学研究中。该方法将扩散函数建模为D(g)=gTDg，其中g是扩散加权梯度方向，D是二阶对称正定张量。这种建模方法虽然简单、稳定，但所采用的二阶对称张量仅能描述体素内的单一平均纤维走向，而无法用于描述存在纤维交叉、分叉等体素内复杂纤维结构。

目前，已有一些方法能够通过更复杂的模型来估计体素内的多个纤维走向，从而解决DTI方法的限制。这些方法主要包括：多张量模型［2-4］、高阶扩散张量［5］、Q-Ball 成像［6-9］、扩散谱成像［10-12］、球面反卷积（Spherical Deconvolution,SD）［13-17］、独立分量分析［18］、混合扩散成像［19］等。其中，基于SD的方法由于具备不需指定纤维分布的数量、计算效率高且能够在低角分辨率成像条件下估计体素内纤维走向分布等优点，得到较广泛的应用，但是基于SD 的方法对噪声非常敏感。为了解决SD 方法对图像噪声敏感的问题，通常有两种改进方案：一是在SD 模型中引入正则化技术；二是在使用SD方法进行体素内纤维结构估计前，先采用有效的降噪方法对DWI 图像进行预处理。非局部均值（Non-Local Mean, NLM）作为一种经典的图像平滑方法已被应用于DWI 降噪［20］，在使用NLM方法对DWI降噪时，只是简单地对一组不同扩散加权梯度方向的DWI 图像分别进行NLM 平滑，往往难以得到较好的结果。

本研究针对NLM 方法直接用于DWI 降噪存在的不足，提出一种面向DWI 的NLM 方法，该方法在对DWI平滑过程中既考虑到不同位置体素间的关联性，同时也考虑到不同扩散加权梯度方向数据间的关联性。原始DWI 数据经本研究提出的方法处理后，再采用SD 方法进行体素内纤维结构估计，可有效提高体素内纤维结构估计结果的准确性。

1 方法

1.1 基于SD的体素内纤维结构估计

SD方法［13］是目前应用较广的体素内纤维结构估计方法，该方法将磁共振扩散信号建模为：

其中，S为扩散加权信号，R为响应函数，F为纤维走向分布函数，⊗为卷积操作。通过构造合适的响应函数R，就可以通过反卷积方法求得纤维走向分布函数F。在该模型下，采用纤维走向分布来表征体素内纤维结构。更多关于SD 的细节信息参见文献［13］。

1.2 经典的非局部均值算法

非局部均值算法由Buades 等于2005年首次提出，该算法基于图像信息的自然冗余进行加权平均降噪。在NLM算法的理论公式中，体素xi的非局部均值NLM(xi)是图像I上所有体素的加权平均，即：

其中，v(xj)为体素xj的灰度值，ω(xi,xj)为体素xj对体 素xi的权重系数 ，ω(xi,xj)∈[0,1] 并且∑xj∈Iω(xi,xj)=1。

在NLM算法的实际应用过程中，对与xi相似体素的搜索通常限定在其邻域Vi中，因此式（2）可写为：

对于在Vi中的每个体素xi，权重ω(xi,xj) 定义为：

其中，Ni和Nj分别为体素xi和xj的邻域，Z(i)=∑jω(xi,xj)为归一化常数，为图像噪声标准差的估计值，h为滤波器参数，距离d定义为：

其中，N为邻域Ni中体素的数量，其值等于邻域Nj中体素的数量，yk和zk分别为邻域Ni和Nj中的第k个体素。对于灰度图像，

1.3 面向磁共振DWI的NLM算法

对于一组DWI 图像，其灰度值（扩散加权信号值）可表示为S(p,g)，p为体素的位置，g为扩散加权梯度方向。考虑角度信息的DWI 的NLM 平滑算法可表示为：

其中，pk为pi的邻域Vpi中体素的位置，G为所有扩散加权梯度方向的集合，权重ωDWI定义为：

其中，Npi和Npk分别为体素pi和pk的邻域，ZD为归一化常数，σ̂D为DWI 图像噪声标准差的估计值，hD为滤波器参数，距离Dis 定义为：

其中，ND为邻域Npi中体素的数量，其值等于邻域Npk中体素的数量，qm和rm分别为邻域Npi和Npk中第m个体素的位置，Δ(S(qm,gj),S(rm,gl))=‖qm-rm‖2lg2(|gjTgl|)。

2 实验与分析

本研究提出的基于NLM平滑的体素内纤维结构估计方法分别在数值仿真数据、仿真实体数据上进行了实验测试，并与不对DWI 图像进行降噪直接采用SD方法［13］以及先利用经典的NLM方法［20］对DWI图像进行降噪后再使用SD 方法（NLM+SD）进行了对比。在数值仿真数据上，通过定量比较不同方法得到的体素内纤维结构的平均角度误差来评价不同方法的准确性，本研究中平均角度误差定义为估计得到的纤维走向与实际纤维走向之间夹角的平均值；在仿真实体数据上，根据已知的纤维路径的走向来推断体素内纤维结构的实际情况，并将不同方法得到的结果与之比较，从而评价不同方法的准确性。

2.1 数值仿真数据实验

实验采用多张量模型［2］的简单形式二张量模型生成一组数值仿真数据，二张量模型是将一个体素内沿扩散加权梯度方向g的扩散加权信号S(g) 建模为两项扩散张量项的和的形式：

其中，b为扩散敏感因子，f为体积分数，D1和D2为两扩散张量，扩散张量可以根据指定的特征向量和特征值进行构造。在数值仿真中，为了简化操作，通常设定S0=1，f=0.5，并假定扩散张量的第二和第三特征值相同，同时令D1和D2的特征值相同。因此只要给定b值，一定数量的扩散加权梯度方向，两个扩散张量的主特征方向，以及张量的特征值，即可生成仿真的扩散加权信号。在本仿真中，给定b值为1 000 s/mm2，30个扩散梯度方向，两个扩散张量的主特征方向如图1a所示，张量的特征值为（1.2,0.3,0.3）×10-3mm2/s。由于数值仿真得到的是无噪声数据，为了测试算法的抗噪性能，向数值仿真数据中加入不同强度的莱斯噪声，从而得到不同信噪比的数值仿真数据。

图1b～图1d给出不同方法在一组信噪比为25的数值仿真数据上的估计结果。从图1b可以看出直接采用SD 方法估计得到的体素内纤维结构呈现出较为混乱的分布，这说明SD 方法极易受到噪声的影响，在噪声强度不高的情况下，也会导致估计结果的误差大大增加。从图1c 和图1d 的结果可以看出，在进行体素内纤维结构估计前对DWI数据进行平滑能够有效提升纤维结构估计结果的质量，减小估计误差。通过对比图1c 和图1d 的结果可以看出，NLM+SD方法一些边缘区域存在模糊现在，例如图1c中红色方框所标注的体素，对于一些实际只有单一纤维走向的体素，NLM+SD方法却得到了两个纤维走向；而对于一些实际上有两个体积分数相同的纤维走向的体素，NLM+SD 方法得到的结果却显示其中一个纤维走向的体积分数较小。而本研究方法在这组仿真数据上并没有显示出这些问题，可以看出本研究方法明显优于NLM+SD方法。

图1 不同方法估计得到的体素内纤维结构Fig.1 Intravoxel fiber structures obtained by different methods

为了定量比较不同方法的估计结果，通过计算图1b～图1d所示估计结果的平均角度误差，得到SD方法的平均角度误差为14.8°，NLM+SD方法的平均角度误差为7.5°，而本研究方法的平均角度误差为4.8°，显然本研究方法得到的结果误差低于其它方法。为了比较不同噪声强度对3种方法（SD方法、NLM+SD方法和本研究方法）的估计结果的影响，在5种信噪比条件下，不同方法估计结果的平均角度误差如图2所示。与SD方法和NLM+SD方法相比，本研究方法的估计结果的平均角度误差约降低25%～30%。

2.2 仿真实体数据实验

实验采用的仿真实体数据来源于MICCAI 2009国际会议，该数据最初是用于人脑白质纤维追踪算法竞赛，关于该组数据的更多细节信息参见文献［21］。该比赛的数据集中包含多组不同参数的DWI数据，本实验中采用其中一组b值为1 500 s/mm2、64个扩散加权梯度方向的数据。

图2 不同信噪比条件下3种方法的平均角度误差Fig.2 Mean angular errors of 3 methods at different signal-tonoise ratios

为了比较经典NLM方法与本研究提出的面向DWI数据的NLM方法，图3a～图3c分别给出未经降噪处理的DWI数据的各向异性分数（Fractional Anisotropy,FA）图、采用经典NLM方法处理后的DWI数据对应的FA图以及采用本研究提出的面向DWI数据的NLM方法处理后的DWI数据对应的FA图。从图中可以明显看出，采用经典NLM方法或本研究方法处理后，相应的FA图的噪声明显降低。与经典NLM方法的结果（图3b）相比，本研究方法的结果（图3c）在存在纤维束的区域内的均匀性更好，在纤维束区域与背景边界更清晰，没有明显的模糊现象。

图3 不同方法得到的FA图和体素内纤维结构Fig.3 Fractional anisotropy maps and intravoxel fiber structures obtained by different methods

为了更好地比较不同方法的效果，选取如图3a～图3c 中红色方框所标识的感兴趣区域，该区域包含两束纤维交叉的情况，分别采用3 种方法（SD 方法、NLM+SD方法和本研究方法）对该区域内的体素进行纤维结构估计，得到的估计结果如图3d～图3f 所示。可以看出，在纤维交叉处SD方法的结果比较混乱，说明SD 方法容易受噪声影响，估计结果的误差较大，而NLM+SD 方法和本研究方法的结果均具有较好的局部一致性。对比图3d～图3f 中红色椭圆标记的区域，该区域的体素内纤维结构应为单一纤维走向，而NLM+SD 方法却得到存在体素内纤维交叉的结果，这说明NLM+SD 方法在存在纤维交叉的区域和不存在纤维交叉的区域的边界处出现了一定的模糊现象，而本研究方法的结果则较少存在这种情况。

3 结论

本研究提出一种基于NLM平滑的体素内纤维结构估计方法，该方法通过建立面向DWI数据的NLM平滑方法，在考虑DWI 数据不同位置的体素间的关联性同时也考虑到不同扩散加权梯度方向的数据间的关联性，对DWI 数据进行NLM 平滑后再采用SD方法估计体素内纤维结构。数值仿真数据和仿真实体数据实验结果表明，与SD 和NLM+SD 方法相比，采用本研究方法得到的体素内纤维结构的平均角度误差更小，且较少存在边缘模糊现象。