基于残差局部幂的分数阶各项异性扩散低剂量CT降噪算法

焦枫媛，桂志国，刘 祎，韩 意

(1. 中北大学 生物医学成像与影像大数据山西省重点实验室，山西 太原 030051；2. 中北大学 信息与通信工程学院，山西 太原 030051；3. 中国人民解放军63961部队，北京 100012)

X射线计算机断层成像(X-ray Computed Tomography, CT)技术在医学领域获得了广泛应用，在CT扫描过程中不可避免地会产生大量辐射剂量，对人体造成伤害[1]. 于是，低剂量CT技术应运而生，如何改善由于剂量降低带来的图像质量下降问题也随之成为当前研究的一个热点，目前的处理方法主要有3类： 投影域降噪[2,3]、统计迭代重建算法[4,5]及后处理算法[6,7]. 其中，后处理算法因其不依赖投影数据且易于操作和推广等优势吸引了国内外学者的关注.

国内外众多学者将经典的图像处理方法改进后应用于低剂量CT图像处理，取得了较好的成果. Chen等[8]在肩部和腹部低剂量CT处理中，应用了大尺度邻域的加权平均算法，有效抑制了噪声和伪影； Kang等[9]在低剂量CT冠脉数据处理中应用了自适应的三维块匹配算法，既有效抑制了噪声，又改善了左心室功能的评估； Chen等[10]提出的三维块匹配算法对于提高图像的对比噪声比效果显著； Wu等[11]提出的级联卷积神经网络算法获得了更高质量的低剂量CT处理图； Wolterink等[12]将生成式对抗网络应用于低剂量CT图像去噪； Shan等[13]在低剂量CT图像处理中将二维训练网络的迁移学习与基于合同路径的三维卷积编码-解码网络相结合，有效保持了图像的精细结构.

基于偏微分方程(Partial Differential Equation, PDE)的方法是图像处理中一类重要方法. 1990年，Perona等[14]提出了著名的各项异性扩散模型，即Perona-Malik (PM)模型，比较以往的各项同性扩散模型，PM模型能在滤除噪声的同时，更好地保持图像边缘. Rudin等[15]于1992年提出了全变分(Total Variatuon, TV)模型，在保持边缘方面有更好的表现. 为降低二阶PDE所产生的块状效应，You等[16]提出了一类四阶PDE，这类算法处理结果图中常常含有斑点. 为克服以上PDE模型带来的问题，Bai等[17]提出了分数阶PM(Fractional-order PM, FPM)模型； Zhang等[18]提出了分数阶TV(Fractional-order TV, FTV)模型； Wang等[19]将FPM与FTV模型加权组合并加以改进提出了分数阶PMTV模型(Fractional-order PMTV, FPMTV)，这类分数阶PDE模型既能克服阶梯效应、斑点效应，又能在一定程度上保持图像的边缘和细节. Hossein等[20]将残差局部幂应用于图像处理，通过残差图像的利用，有效提高了PDE模型降噪过程中细节和纹理的保持.

由于低剂量CT图像中存在大量组织结构和细节，在去除噪声的同时，最大限度地保持图像细节对于医疗诊断具有重要意义，受文献[20]启发，本文将局部残差幂作为纹理检测算子，引入FPM模型的扩散系数中，使其与梯度一起自适应调节模型扩散系数，使得获得的低剂量CT处理结果图在去除噪声的同时更好地保持了图像纹理和细节信息. 通过头部体模和实际临床数据实验验证了本文算法的有效性.

1 分数阶微积分

Grümwald-Letniklv (G-L)和Riemann- Liouville (R-L) 定义是分数阶微积分各类定义中应用最广泛的，本文将G-L定义应用于低剂量CT图像处理中，信号f(x)的α阶微分的G-L定义如下[18]

(1)

2 基于残差局部幂的分数阶各项异性扩散算法

2.1 整数阶各项异性扩散算法

整数阶各项异性扩散算法，即PM模型的表达式为

(2)

式中，扩散系数函数c可表示为

c(|∇f|)=1/[1+|∇f|2/k2] ,

(3)

式中：k为梯度阈值. PM模型虽然取得了很好的去噪效果，但是降噪后的图像经常会出现块状效应.

2.2 基于残差局部幂的分数阶各项异性扩散算法

为了克服阶梯效应，Bai等[17]在2007年提出了分数阶各项异性扩散方程，即FPM模型，其能量函数可表示为

(4)

式中：Ω为图像支撑，且p(|Dαf|)≥0是一个增函数，其与扩散系数关系可表示为

(5)

为了解决上述问题，更好地辅助临床诊断，本文将纹理检测算子添加到分数阶各项异性扩散方程的扩散系数中. 由文献[20]可知，残差局部幂可作为一种纹理检测算子，首先定义残差图像

fR=f0-f=fs+fn,

(6)

式中：fs表示纹理和细节成分；fn表示噪声成分，由于两者具有不相关性，因此

PR=Ps+Pn,

(7)

式中：PR、Ps和Pn分别为fR、fs和fn的局部幂，将PR做如下归一化处理后可得

(8)

由此，分数阶各项异性扩散模型中扩散系数可写为

(9)

式中：k1为常数.

根据文献[17]中分数阶微分的离散化方法，将(4)式中的分数阶微分替换为卷积积分后可表示为

(10)

式中：

根据梯度下降法，可以得到

(11)

将(vα*f)x和(vα*f)y离散化获得

(12)

记

式中：ε为很小的整数. 基于局部残差幂的分数阶各项异性扩散方程的离散化形式可表示为

[Dx(i+k,j)+Dy(i,j+k)].

(14)

2.3 算法描述

1) 通过FPM模型将噪声和相关纹理进行分离. 将低剂量CT图像fs与FPM处理结果图f0的差作为残差图，当所得残差的方差等于低剂量CT图像噪声方差的β倍时，FPM算法停止迭代.

2) 计算残差局部幂. 将步骤1)停止迭代后获得的残差图通过高斯滤波滤除部分噪声后，再应用式(8)获得归一化的残差局部幂.

3) 将基于残差局部幂的FPM模型应用于低剂量CT图像. 将步骤2)获得的残差局部幂加入分数阶各项异性扩散模型扩散系数中，根据式(9)～式(14)进行计算，获得最终处理结果图.

3 实验结果及分析

为了验证所提算法的有效性，本文对Shepp-Logan头部体模和实际临床数据进行仿真实验. 头部体模大小为256 mm×256 mm，实际临床数据为DICOM医学图像库中大小为512 mm×512 mm 的腹部高剂量CT与低剂量CT图. 本文与FPM算法[17]、改进的非局部均值(improved Non-Local Means, INLM)算法[21]和FPMTV算法[19]进行了比较. 实验环境如下： 操作系统是64位Windows 10，处理器是Intel Core(TM)i7-8565U CPU @ 1.80 GHz 1.99 GHz，内存为32 GB. 比较算法的参数根据相应参考文献选择，本文算法参数设置为β=1.4，α=1.03，Δt=0.25，迭代次数在头部体模实验中设为2 200次，在腹部数据中设为55次.

图 1 头部体模各算法处理结果图及局部 感兴趣区域放大图对比Fig.1 Comparison of the denoised images and their zoom-in ROIs processed by a variety of algorithms on Shepp-Logan phantom

图 1 给出了各种算法应用于头部体模所获得的处理结果图及局部感兴趣区域放大图对比. 从图 1 可以看出，FPM和FPMTV算法处理结果图不同程度地模糊了图像的边缘和细节，INLM算法处理结果图在噪声较多的区域仍然有明显的伪影留存，本文算法在滤除噪声的同时，很好地保持了图像的边缘和细节.

图 2 给出了各种算法处理腹部低剂量CT图的结果对比，在实际临床数据中，高、低剂量图中均含有噪声点，但低剂量图中含有更多斑点噪声，不利于腹部组织结构的显示. 通过比较各算法处理结果图可以看出，本文算法在去除噪声点的同时更好地保持了基本的组织结构.

图 2 腹部低剂量CT图各算法处理结果对比Fig.2 Comparison of the denoised images processed bya variety of algorithms onabdominal low-dose CT image

图 3 给出了腹部低剂量CT图与不同算法处理结果图的差值图. 从图3(a)和图3(c)可以看出，FPM和FPMTV算法在噪声被滤除的同时，部分边缘和结构也被滤除，从图3(b)和图3(d)可以看出，INLM和本文算法在保持边缘和细节方面明显优于前两个算法，进行细微比较则可以发现本文算法残差图中几乎看不到残留的边缘和细节信息.

为了对本文算法的有效性进行定量验证，在Shepp-Logan头部体模实验中，采用峰值信噪比(Peak Signal to Noise Ratio，PSNR)和结构相似度(Structural Similarity，SSIM)对各算法进行定量描述. 在实际临床数据实验中，对于病人呼吸和运动等因素带来的高、低剂量CT图难以完全匹配问题，采用标准差(Standard Deviation, STD)进行客观评价. 各评价指标定义为

(7)

(8)

式中：

σuoriginalu=Cov{uoriginal,u}=

(9)

图 3 腹部低剂量CT图与不同算法处理结果图的差值图对比Fig.3 Comparison of the image subtractions between the abdominal low-dose CT image and the denoised images by a variety of algorithms on abdominal low-dose CT image

表 1 为FPM、INLM、FPMTV以及本文算法分别应用于头部体模而获得的处理结果图相对于原始图像的PSNR和SSIM值的比较. 由表1很容易看出，本文算法的PSNR值最高，说明本文算法处理结果图中含有最少的噪声； 本文算法的SSIM值最高，说明本文算法处理结果图与原图最接近.

表 1 头部体模各种算法的客观评价Tab.1 The objective evaluation of a variety of algorithms on Shepp-Logan phantom

图 4 为低、高剂量腹部CT图像以及不同算法处理结果图的4个ROIs的STD值柱状图，图 4 左上角的腹部图用方框标注了4个ROIs. 通过柱状图可以直观地看出，在4个ROIs中，本文提出的FPMLP算法对应的STD值均最接近高剂量CT图，这表明本文算法处理结果图与高剂量CT图最接近. 结合图 1～图 4 和表 1 可知，在视觉效果和定量评价方面，本文算法都优于其他算法，既有效降低了低剂量CT图像的噪声和伪影，又较好地保持了CT图像的纹理、细节等结构特征.

图 4 低、高剂量CT图像以及各算法处理结果图在腹部图方框标注的4个ROIs的STD值柱状图

4 结 论

本文提出了一种基于残差局部幂的各项异性扩散低剂量CT降噪算法. 该方法将残差局部幂作为纹理检测算子加入分数阶各项异性扩散模型的扩散系数中，充分利用了残差图像中的纹理、细节和边缘信息，有效增强了分数阶各项异性扩散算法的降噪性能. 基于Shepp-Logan和实际腹部医疗数据的仿真实验结果表明，本文算法在抑制噪声和保持图像纹理、细节和边缘方面均优于传统的低剂量CT降噪算法.