基于结构张量和各向异性平滑的DTI去噪

刘帅奇,李鹏飞,安彦玲,扈 琪,赵 杰

1(河北大学 电子信息工程学院,河北 保定 071000)2(河北省数字医疗工程重点实验室,河北 保定 071000)3(河北省机器视觉工程技术研究中心,河北 保定 071000)

1 引 言

扩散张量成像(DTI)已经演变成为具有非侵入性的一种表征活体组织结构的主要技术[1].这种表征技术是基于所测量的扩散张量的特征参数,其中包括通过特征值的差异来评估组织结构完整性,并且通过特征向量的方向来评估结构特征.然而特征参数对噪声极其敏感,因此随机噪声会影响结构特征的系统偏差评估.由于DTI数据信噪比较低,且噪声不仅会引起信号的随机波动,还会导致信号与数据的偏差,特别是进行纤维结构分布的分析时,被噪声污染的张量通常表现为方向排列上杂乱、不规则,使得实验得到的纤维结构不够平滑,甚至出现错误的结果,大大的限制了DTI的应用.因此抑制扩散张量图像中产生的噪声有利于减少这些系统偏差,从而提高DTI组织表征的准确性.

对于所有矢量或张量图像平滑来说,减少扩散张量图像中的噪声需要两个基本要求:恢复包含在张量数据中的方向信息和保留结构边界.为了实现这两个目的,Parker等人应用了Perona和Malik提出的非线性扩散滤波技术[2]来减小扩散张量图像中的噪声.虽然可以有效地降低扩散各向异性的有偏估计,但是平滑对张量方向信息的影响只是停留在定性证明的阶段.Pajevic等人应用样条拟合法从离散DTI数据中生成一个连续、平滑的张量场[3],但是这种方法不能有效地保留数据中的结构边界.Poupon等人提出的基于马尔可夫模型[4]和随机弛缓标记[5]的用于恢复主扩散方向(与最大特征值相对应的方向)的正则化方法,也因其边界保留能力而存在一定的局限性.

结构张量技术可以很好地将结构信息突出的部分和结构信息微弱的部分进行区分,例如区分图像中边缘轮廓等细节与平坦光滑的部分,这样可以大大地方便类似医学图像等对边缘保留要求比较高的图像去噪.Bigun和Granlund等人将结构张量应用于图像方向场的估算[6]等方面,利用结构张量快速提取图像中的结构信息,例如图像中目标边缘、目标形状角点特征等.Weickert和Brox对结构张量的线性高斯滤波技术做了很大的改进,利用基于偏微分方程形式的扩散方程的方法将线性滤波过程变成了非线性[7],这样改进后很大程度上提高了边缘检测和角点检测的精度.

迄今为止,国内外科研人员已经提出了大量的图像平滑技术来减少图像中的噪声.其中包括B-样条拟合[8]以及变分原理[9]、差分几何[10]、变换域去噪算法[11]和经验分解去噪算法[12].然而,这些平滑技术的数据库尚未实现其实用性,部分原因在于数值实现方面有些耗时,特别是计算复杂度会随着加权方向数量的扩大而增加,或者是对其实用价值而言,缺乏严格的体内DTI数据证明.各向异性平滑的概念最初由Weickert提出,用于增强标量图像中的流状结构[7].各向异性平滑与传统滤波方法的不同之处在于它的各向异性,即在不同局部区域、不同方向上,滤波与平滑的强度存在一定差异.它与图像的局部结构(灰度梯度)相关,在灰度变化较大的区域或方向上,由于灰度梯度较大、扩散作用较小,从而最大限度地保留这些局部细节特征;若灰度变化较小或只有孤立噪声点的区域,则采用较强的平滑处理以实现抑制噪声、保留边缘与细节的目的.这一概念的提出使得平滑技术有了质的飞跃,其对减少图像中的噪声非常有效,并且使得结构边界更加平滑.

在本文中,我们结合结构张量和各向异性平滑技术提出了一种新型DTI去噪方案,首先利用结构张量区分DTI图像中的边缘结构信息和平坦光滑信息,然后对图像结构内部的平坦光滑区域进行各向同性平滑,对结构边界区域进行各向异性平滑,最终得到去噪后的DTI图像.

2 基本理论

2.1 结构张量

结构张量是从图像像素中提取结构信息的张量,通常应用于图像和视频处理[13,14].定义Ω⊂R2表示一幅二维图像数据矩阵,该二维图像的灰度I:Ω→R,那么这幅图像中的每一个像素的初始结构张量就是一个2×2的对称矩阵场,通过该矩阵场可以求出包括边缘、方向等有关灰度图像的信息.

初始结构张量ST(i,j)定义为

(1)

通常,我们在平滑处理时将结构张量定义为

(2)

其中Gρ是标准差为ρ的高斯核函数.因为和高斯核函数的卷积是一个线性过程,很多文献都称这种经典的结构张量为线性结构张量[15-17].因为抑制了噪声水平,经过高斯核函数处理过的线性结构张量能够更好的体现出方向信息和边缘信息.

对结构张量的每个像素信息即对公式(2)进行特征分解,可以得到两个特征值λ1≥λ2,其对应的特征向量为v1和v2,v1近似平行于最小梯度的方向,而v2与该方向正交并且大致平行于最大梯度的方向,每个特征值的大小是其对应的特征向量方向上的梯度量.

2.2 各向同性平滑

高斯滤波是各向同性平滑的经典模型.其根据所选择的高斯核函数来控制向周围扩散的比重参数[10].高斯滤波对于去除高斯白噪声具有良好的性能.均值为零、方差为σ的一维高斯函数可以表示为:

(3)

由于DTI图像一般是二维或三维的数据,所以通常采用二维高斯函数作为各向同性滤波器,其表达式为

(4)

其中A为归一化参数,σ为噪声方差,其决定了滤波器的去噪效果,方差越大,去噪效果越好.但是这种滤波方法也存在一定的不足之处,那就是在去除图像噪声的同时产生图像细节和边缘模糊现象.

2.3 各向异性平滑

为了增强标量图像中的流状结构,Weickert提出了使用偏微分方程平滑图像的各向异性扩散滤波技术[7]

(5)

(6)

其中⊗代表外积运算,参数ρ是高斯核的标准差,决定了梯度张量的空间尺度.G的特征值(λ1,λ2,λ3,以降序排列)和特征向量对应于图像强度梯度在空间尺度ρ上的大小和方向.与最大特征值相关联的特征向量是最大强度梯度的方向,与最小特征值相关联的特征向量是最小强度梯度的方向.通过将特征值不成比例地改变尺度(但保留特征向量),使得对应于最小强度梯度方向的特征值是最大的,以这种方式构造的结构张量可以有效地提升图像中沿边缘切线方向区域的平滑程度.

3 基于结构张量和各向异性平滑的DTI去噪

在本节中,详细说明了本文提出的基于结构张量和各向异性平滑的DTI去噪算法的实现过程.算法的主要步骤包括:(1)结构张量对含噪DTI图像的提取;(2)对含噪DTI图像进行各向异性平滑.具体来说,首先利用结构张量从DTI图像像素中提取结构信息,并得到图像的均匀平坦区域和边缘轮廓区域.然后在结构内部即平坦区域各向同性地平滑图像,在结构边界即边缘轮廓区域各向异性地平滑图像.最终得到去噪后的DTI图像.

3.1 结构张量对含噪DTI图像的提取

我们将2.1节中的结构张量矩阵E记作

(7)

其中Ix,Iy就是对原图像在x和y方向求得的偏导.然后我们求出矩阵E的行列式K=det(T)和迹H=tr(T),并且根据行列式K和迹H的关系来区分图像的均匀平坦区域和边缘轮廓区域.区分原则为:当H=0时,认为结构张量提取的是图像的均匀平坦区域:当H>0 &K=0时,则认为是边缘轮廓区域.

3.2 对含噪DTI图像进行各向异性平滑

(8)

为了在结构内部各向同性地平滑图像,在结构边界各向异性地平滑图像,结构张量T由上述梯度张量G构成,将T定义为G的归一化倒数.vg1,vg2,vg3表示G的特征向量,λg1,λg2,λg3表示G的相应的特征值,T的特征向量和特征值是:

vti=vgi,i=1,2,3

(9)

λti=1/λgi,i=1,2,3

(10)

根据上述公式可以看出,在均匀区域(内部束状结构)中,结构张量与梯度张量的三个特征值的大小是相同的,并且沿着所有方向的各向同性平滑结果是相似的.在结构边界,垂直于边界方向(较大强度梯度的vt1方向)上的结构张量的特征值是较小的,沿着边界的方向(较小强度梯度的vt3方向)上的结构张量的特征值是较大的,因此图像结构在沿着结构边界的切线方向上比垂直方向更平滑(各向异性平滑).为了实现在均匀区域和结构边界上产生同样的平滑,结构张量的特征值被归一化为常数C,即λt1+λt2+λt3=C,因此无论是各向同性还是各向异性,平滑的总量在整个图像上是相同的.

图1 算法设计流程图Fig.1 Flow chart of our algorithm

为了总结本文所提出的算法的实现过程,对算法步骤给出了简明地描述,算法设计流程图如图1所示.

基于结构张量和各向异性平滑的DTI去噪算法步骤:

步骤1.构造结构张量矩阵E,求出矩阵E的行列式K和迹H,然后根据行列式K和迹H的关系来区分图像的均匀平坦区域和边缘轮廓区域.

步骤3.使用公式(9)和公式(10)构造结构张量T.结构张量T应为一个3×3的对称正定矩阵,然后分析并计算它的特征值和特征向量,从而实现在结构内部的均匀平坦区域各向同性地平滑图像,在结构边界的边缘轮廓区域各向异性地平滑图像,最终得到去噪后的DTI图像.

4 实验结果和分析

本文算法及对比算法使用MATLAB编写,实验平台为Intel Core i5 3.20GHz CPU和4.00 GB RAM的工作站.本文首先使用真实大脑的数据来验证结构张量是否能很好的提取图像的平坦区域和边缘区域.在实验中,先求出第3节中提出的结构张量矩阵T的行列式K和迹H,并根据行列式K和迹H的区分原则得出平坦区域和边缘区域的分割结果.图2给出了利用结构张量区分出来的平坦区域和边缘区域.这表明结构张量可以有效提取图像中像素并加以区分.

利用结构张量将图像区分并且得到平坦区域和边缘区域之后,对含噪DTI图像应用本文所提出的各向异性平滑算法进行处理.为了验证本文所提出的各向异性平滑算法在边缘保留和去噪方面的良好性能,将本文所提的各向异性平滑算法与文献[17]中基于复剪切波及复扩散各向异性滤波(CST-CDAF)的DTI去噪算法和文献[18]中基于半隐式(S-I)方案的各向异性平滑算法进行比较.

图2 结构张量提取结果Fig.2 Extraction results of structure tensor

首先用模拟的DTI数据集来测试本文所提出的各向异性平滑算法的有效性.数据集由两块不同的“纤维”束组成,每个块包含沿着不同取向的“纤维”片段.“纤维”的平均扩散系数(D)为0.7×10-5cm2/s,部分各向异性(FA)为0.9,与人脑体内测量值相似.

上述各去噪算法的去噪结果如图3所示.其中图3(a)为模拟 “纤维”分布结构的DTI图像,图3(b)为对模拟DTI数据添加零均值、SD=0.10的高斯噪声的结果,图3(c)为CST-CDAF算法的去噪结果,图3(d)为S-I算法的去噪结果,图3(e)为本文所提算法的去噪结果.可以看出在平滑后,由于噪声导致的主扩散方向(PDD)的变化程度明显降低,并且在任何“纤维”边界附近没有观察到模糊现象,这也证明了我们所提出的算法在边缘保留和去噪方面的有效性.

图3 各向异性平滑算法对模拟DTI数据的影响Fig.3 Influence of anisotropic smoothing algorithm onsimulated DTI data

为了测试本文所提出的各向异性平滑算法对人脑DTI数据的去噪性能,使用磁共振扫描仪从健康的受试者大脑中获得扩散加权图像.DTI数据是基于自旋回波单次激发EPI成像获得的,每次扫描中的成像参数b=1000sec/mm2,并且扫描体积为250×250×120mm2.扫描所获取的数据矩阵的大小为128×128×34,然后运用插值法将其转换为256×256×34.在扫描过程中,本文进行多次重复扫描并求其平均值,最终产生信噪比为SNR≥70%的数据集.选择数据集中的两个代表性区域分别应用上述各去噪算法,并且与模拟数据中的参数保持一致.两个代表性区域及其PDD结构如图4所示,第一个区域含有白质,其中存在具有不同方向的纤维束.第二个区域包含胼胝体的一部分,其中可以看见大部分平行的神经纤维.

图4 数据集中两个代表性区域及其PDD结构.左为扫描后获得的人脑扩散加权图像,右上方为人脑中含有白质区域的PDD结构,右下方为人脑中含有胼胝体区域的PDD结构Fig.4 Two representative regions of the dataset and their PDD structure.The left image is a diffusion weighted image of the human brain obtained after scanning,the top right image is the PDD structure containing the white matter region in the human brain,and the bottom right image is the PDD structure containingthe corpus callosum region in the human brain

上述各去噪算法对人脑中两个代表性区域的去噪结果如图5和图6所示.从图中的实验结果可以看出,本文所提出的各向异性平滑算法明显降低了噪声对两个区域中主扩散方向(PDD)的影响,并且不会模糊结构边界.

图5 各向异性平滑算法对人脑中含有白质区域的DTI数据的影响Fig.5 Influence of anisotropic smoothing algorithm on DTI data containing white matter regions in human brain

为更加准确的说明本文所提算法的有效性,我们将上述各向异性平滑算法对PDD的影响进行定量评估.平滑后PDD的改善程度通过平均角度差改变的百分比来衡量,其定义为

(11)

其中θ0,θN分别是平滑之前和N次平滑之后DTI数据中PDD的角度差.上述各去噪算法对人脑中两个代表性区域PDD的平均角度差变化的影响如图7所示,图中显示了在大脑两个代表性区域中PDD的平均角度差随迭代次数变化的情况.从图中可以很容易看出基于CST-CDAF去噪算法出现PDD的平均角度差随着迭代次数的增加而增大的情况,这表明在平滑过程中出现了结构模糊现象.基于S-I去噪算法和本文所提算法均呈现出PDD的平均角度差随迭代次数的增加而减小的趋势.对本文所提的各向异性平滑算法的定量分析表明,两个代表性区域PDD的平均角度差均随着迭代次数的增加而减小,其中迭代次数N= 1、2、4和 8对应的两个区域的平均角度差改善的百分比分别是27%、40%、60%和92%.综上所述,定量评估的结果表明本文所提的各向异性平滑算法在去噪和保留结构边界等方面优于CST-CDAF算法,并且可以得到和S-I算法不相上下的结果.

图6 各向异性平滑算法对人脑中含有胼胝体区域的DTI数据的影响Fig.6 Influence of anisotropic smoothing algorithm on DTI data containing corpus callosum area in human brain

图7 定量比较三种不同的平滑算法对人脑中两个代表性区域的实验结果Fig.7 Experiment results of the quantitative comparison between three different smoothing algorithms on the two representative regions of the human brain

5 总 结

在本文中,我们在结合结构张量和各向异性平滑技术的基础上提出了一种新的DTI去噪算法.该算法将结构张量区分图像结构的特性和各向异性平滑的边缘保留能力完美的结合在一起.基于模拟和体内数据的实验表明,本文所提出的去噪算法在边缘保留和去除噪声方面都显示出了很好的性能,对医学图像后期处理和病理分析产生深远的影响.