0 引言
为了解决CT检查导致的辐射剂量过高问题，低剂量CT扫描逐渐成为临床应用热点.然而，低剂量CT扫描会导致图像质量退化，特别是LDCT图像中的噪声会严重影响诊断精度.因此，去除噪声是改善低剂量CT图像质量的重要任务.低剂量CT扫描所引起的噪声主要是X射线量子噪声[1].针对图像中的细节和噪声同属高频成分难以分离的经典难题，小波变换等多尺度分析方法呈现出一定的优势[2].近几年，许多学者已提出大量的小波域去除量子噪声方法[3-4].主要包括两种思路：其一，假设其服从Gaussian分布，利用阈值法去噪，如小波硬阈值、软阈值和贝叶斯阈值法等[3].其二，假设其近似服从Poisson分布.首先利用方差稳定变换将其转化为近似服从Gaussian分布的噪声，再利用去除Gaussian噪声的常用方法进行处理[4-7].随着多尺度分析方法的发展，一些学者提出了基于剪切波变换的图像去噪算法.例如，Dan等[8]提出的基于剪切波变换的贝叶斯最大后验估计方法，能有效去除因小波去噪产生的伪吉布斯效应；胡海智等[9]提出的全变差正则化剪切波域收缩方法，能更有效地去除Gaussian噪声并保留图像中的细节信息等等.为了去除低剂量CT图像中的量子噪声，笔者假设量子噪声近似服从Poisson分布，并利用Anscombe变换，将原图像中的噪声转化为一个具有近似常数方差的Gaussian噪声.然而，由于传统Anscombe变换对低信噪比的图像其噪声方差无法趋于稳定[10]，致使去噪效果欠佳.为了解决这一问题，文献[11]通过最优化Anscombe逆变换以改进去噪效果；Zhang等[12]通过改进Anscombe变换，使当Poisson噪声参数λ>0.1时，其噪声方差可以趋于稳定.但是，当λ≤0.1时，该方法仍不能使噪声方差稳定到一个常数值.基于前期研究[7]可知，当Poisson噪声参数λ≤0.1时，如在应用Anscombe变换的基础上，再利用阈值法去除量子噪声可以取得更好的效果，但其去噪效果依赖于噪声方差值的估计精度.
鉴于此，笔者针对低信噪比的剪切波高频系数子带，提出基于剩余自相关功率的噪声方差估计方法，并结合贝叶斯最大后验估计提取出非噪声高频系数，从而有效地提高低信噪比图像的去噪效果.
1 剪切波变换
虽然小波变换作为多尺度分析的代表方法被广泛应用于图像去噪，但由于理论上小波对具有线状或面状奇异的高维目标函数而言，并不能达到最优稀疏表示.为了克服小波变换的局限性, 基于复合伸缩小波理论提出的剪切波变换应运而生.由于非下采样剪切波变换(non-subsampled shearlet transform, NSST) 不仅具有较好的稀疏表示特性，而且利用NSST对图像进行分解，能更好地保留原图像中的细节信息.因此，笔者选用非下采样离散剪切波变换，其实现流程如图1所示.首先进行多尺度分解.用拉普拉斯金字塔分解算法把尺度为j-1的近似系数
分解为尺度为j的低频系数
和高频系数
非下采样剪切波变换在每次多尺度分解前，在分解滤波器上采样，使每层的子图像大小一致.多尺度分解由小波基实现，通过尺度分解从而得到不同尺度的子图像.然后，对每一尺度的高频系数
进行方向局部化.计算在伪极坐标系上的离散傅里叶变换.然后，基于改变位置的窗函数即带通滤波器,得到各个方向的高频信息.其中，各个尺度的方向数由向量N=[n1,n2,…,nj]决定.剪切波分解后在尺度J(j=1,2,…,J)上有2nj个方向.图2表示一个原图像及在尺度为1时的剪切波分解示例.其中图2(a)表示原图像，图2(b)、(c)分别表示尺度为1时的低频系数和8个方向的高频系数.
图1 剪切波分解算法框架
Fig.1 The block diagram of shearlet decomposition
图2 剪切波分解的一层图像
Fig.2 The shearlet coefficients at one scale
2 改进的噪声方差估计算法
虽然Donoho提出的噪声方差计算方法(MAD方法)被广泛应用,但因其是针对小波域估计高斯噪声方差而设计，存在一定误差.为此，剩余自相关功率(RAP)[10]与贝叶斯最大后验估计相结合的改进噪声方差估计精度的方法，其实现过程如下.
2.1 贝叶斯最大后验估计
设
g=f+w，
(1)
式中：f和g分别表示原图像和加噪后的图像；w表示均值为0的Gaussian白噪声.
设
y=x+n，
(2)
其中，y和x分别表示含噪声和无噪声图像的系数；n是噪声的剪切波系数.剪切波域去噪即从观测系数y中估计无噪声图像系数x的近似值
基于贝叶斯最大后验估计计算该近似值的步骤如下.
(1)用Donoho提出的估计方法计算噪声方差近似值:
(3)
其中，yi,j表示选取的高频子带系数.
(2)计算提取非噪声剪切波高频子带系数的阈值如下：
(4)
(5)
(6)
(3)用贝叶斯最大后验估计法提取高频子带中非噪声剪切波系数近似值：
(7)
2.2 噪声方差估计
对CT图像加入均值为0的Gaussian白噪声，用非下采样剪切波对图像进行尺度为4的分解，每层选定一个方向.由于随着分解层数的增加，高频子带的图像信息逐渐减少，而噪声信息逐渐增多，且第4层高频子带的信噪比最低，故笔者分别用MAD方法和剩余自相关功率方法计算含噪图像第4层高频细节子带的噪声方差，并对实验结果进行比较，其实现步骤如下.
(1)选取噪声方差的候选值.利用式(3)计算所选高频子带的噪声方差σm，并且以该值为中心,间隔为0.1, 选取30个噪声方差候选值∑=[σ1,σ2,…,σm,…,σ30].
(2)计算所选高频子带所需的阈值.针对选定的每一个噪声方差候选值σm，利用式(4)所示的贝叶斯最大后验估计方法计算阈值[8].
(3)计算剩余误差.利用贝叶斯最大后验估计法估计第m个噪声方差对应的去噪图像
并计算其剩余误差：
(8)
(4)计算剩余自相关功率.计算第m个噪声方差候选值对应的剩余误差自相关
然后，计算剩余自相关功率.
(9)
其中,N是自相关中点的数量.当真正的噪声方差为20时，基于不同的噪声方差候选值计算出其对应的Pm，其对应关系如图3所示.从图3可以看出，其对应的Pm在接近真正的噪声方差时有一个突变，且之后Pm趋于稳定.
(5)计算噪声方差.为了捕捉Pm的这种特性，记
Dm=Pm+1-Pm.
(10)
不同的噪声方差候选值与Dm的对应关系如图4所示，从图4可以看出，真正的噪声方差即当Dm达到最大值时对应的σm右侧一个较小的值，这个值的估计公式为:
mmax=argmmaxDm.
(11)
(12)
σ(RAP)=σm*.
(13)
图3 当噪声方差为20时的剩余自相关功率
Fig.3 Residuals autocorrelation power (RAP) when the true variance is 20
图4 对应图3不同的RAP由式(10)所得的Dm
Fig.4 Difference of RAP, Dm in(10),for the RAPs in Fig.3
为了验证该方法的准确性，对常规剂量的CT图像加入均值为0、方差不同的Gaussian白噪声进行一系列的仿真实验，实验结果如表1所示，表示同一幅CT图像加入的Gaussian白噪声方差分别为5、10、15、20、25时的实验结果.
表1 同一图像Gaussian白噪声方差不同时实验结果对比
Tab.1 Comparison of results in the same image with different White Gaussian noise variances
σ(Gaussian)＾σ(本文方法)＾σ(MAD)55 117 44109 979 931514 7912 572020 1715 302525 2018 04
由表1的实验结果可以看出，对不同的Gaussian白噪声方差，本文算法的估计精度均高于文献[3]中的MAD方法.
3 本文算法
笔者提出的自适应低剂量CT图像质量改善算法，其算法框图如图5所示，主要包括5个步骤.
图5 本文算法框图
Fig.5 The diagram of the proposed algorithm
步骤1：Anscombe变换.首先采用Anscombe变换将低剂量CT图像中的量子噪声转变为具有近似常数方差Gaussian噪声[10].变换公式为：
(14)
式中:
表示噪声图像；
表示Anscombe变换后的图像.
步骤2：剪切波分解.将Anscombe变换后的图像进行4层剪切波分解，每层分解为8个方向，得到1个低频图像和32个高频图像.
步骤3：提取非噪声系数.由于随着分解层数的增加，高频子带的信噪比逐渐降低，故利用提出的噪声方差估计方法对信噪比较低的第3、4层高频子带估计其噪声方差.然后基于噪声方差估值，结合贝叶斯最大后验估计方法提取非噪声高频系数.
步骤4：重构图像.基于步骤3获得的非噪声高频系数和低频系数进行剪切波逆变换，获得去噪图像；对去噪后的图像实行Anscombe逆变换，得到重构图像.公式为：
(15)
式中：
表示去噪后图像；
表示重构图像.
4 仿真实验
笔者从定量评价和视觉效果来验证提出算法的有效性.采用PSNR(peak signal to noise ratio)值和MSSIM (mean structure similarity)值[11]为定量评价指标:
(16)
(17)
其中，f表示原图像,
表示去噪后的图像.式(16)中,fj和
分别表示f和的第j个局部窗中的子图像； M为局部窗子图像的个数.式(17)中，
为亮度比较函数；
为对比度比较函数；
为结构比较函数；α、β、γ分别是用来调整亮度、对比度以及结构信息的权重，使用大小为11×11的加权Gaussian窗函数, 并设置α=β=γ=1.
4.1 定量评价
对取自郑州大学第一附属医院的一系列含有病变的、大小为512×512的常规剂量的胸部CT图像加入参数λ≤0.1的Poisson噪声，并进行仿真实验，采用PSNR和MSSIM来定量评价去噪效果,实验结果如表2和图6所示.
表2表示对10幅常规剂量的CT图像添加强度参数λ=0.01的Poisson噪声，利用本文方法与文献[3]方法对噪声图像进行去噪的实验结果.
表2 当Poisson强度λ=0.01时的PSNR和MSSIM比较
Tab.2 Comparison of PSNR and MSSIM when the Poisson noise intensity is 0.01 dB
图像噪声图像文献[3]方法本文方法PSNRMSSIMPSNRMSSIMPSNRMSSIM118 560 66720 030 71529 580 961218 630 66520 180 71532 210 969318 720 67120 150 72429 220 969418 440 67019 990 71029 470 969518 520 66219 920 72330 160 970618 450 66320 560 71230 270 969718 670 67019 160 71630 250 968818 530 66520 420 71430 830 969918 780 64519 500 71431 260 9681018 560 67220 430 71529 320 968
图6表示对一幅CT图像添加强度参数λ为0.005、0.008、0.010、0.050、0.080和0.100的Poisson噪声，本文方法与文献[7]方法对噪声图像去噪的实验结果.图6(a)表示随着噪声量的变化，两种算法的重构图像与噪声图像的PSNR折线图，图6(b)表示随着噪声图像的MSSIM变化，两种算法的重构图像的MSSIM折线图.
图6 不同Poisson噪声强度时的PSNR和MSSIM对比图
Fig.6 Comparison of PSNR and MSSIM in image with different Poisson noise intensity
由表2可知，本文方法与文献[3]方法相比，PSNR和MSSIM分别提高了52.2%和34.9%.从图6可以看出，对参数λ≤0.1的Poisson噪声，本文算法所得重构图像的PSNR和MSSIM均高于文献[3]方法；本文方法在噪声量较多的情况下，去噪效果尤为突出.
表2和图6的实验结果表明，对参数λ≤0.1的Poisson噪声，无论是同一噪声强度下的不同图像还是多种噪声强度下的同幅图像，与文献[3]方法相比，本文算法所得重构图像的PSNR和MSSIM均有所提高.由此得出，本文算法具有较强的自适应性和鲁棒性.
4.2 视觉效果
本文视觉效果评价采用在日本千叶大学先端医、工学研究中心，以体膜为扫描对象，进行低剂量扫描所采集的降低不同剂量的多幅实际CT图像.低剂量扫描的主要原理是利用放射线剂量与管电流的线性关系，当保持其他参数不变时，通过降低管电流和管电压来降低放射剂量.图7(ai)(i=1,2)所表示图像的扫描参数设置分别为管电压120 kV，扫描层厚为0.5 mm，管电流分别为100 mA和150 mA；图7(a3)表示图像的取像参数设置分别为管电压为120 kV，管电流为200 mA，扫描层厚为1 mm.
从图7(ai)(i=1,2,3)可知，随着管电流的减少，噪声量越来越多.分别用本文方法与文献法对这些图像进行仿真实验，实验结果分别由图7(bi)和图7(ci)(i=1,2,3)表示.图7(bi)(i=1,2,3)为利用文献[3]方法所获得的去噪图像，图7 (ci)(i=1, 2,3)分别表示本文方法的重构图像.由图7可以看出，本文算法不仅能有效去除噪声，而且与文献[3]方法相比能更好地保留图像中的边缘信息和纹理信息.
上述实验的定量分析和视觉效果表明，本文算法具有自适应地去除参数λ≤0.1的Poisson噪声的能力，而且该方法有效提升了重构图像的质量.
图7 低剂量CT图像去噪效果对比
Fig.7 The comparison of de-noising effects for low-dose CT
5 结论
笔者利用剪切波良好的稀疏特性，将Anscombe变换与改进的自适应噪声方差估计算法相结合，提出了一种有效改善低剂量CT图像质量的方法.对含有λ≤0.1的Poisson噪声的CT图像进行的定量评价.结果表明，本文方法与文献[3]方法相比，其PSNR平均提高52.2%，MSSIM提高34.9%.同时，利用本文方法和文献[3]方法对多幅实际低剂量CT图像进行去噪的对比实验结果也表明，本文方法较文献[3]方法可以更好地保留图像细节信息.因此，本文算法具有自适应能力强、重构图像质量高的特点，具有一定的临床应用价值.
[1] SPERL J, BEQUE D, CLAUS B, et al. Computer-assisted scan protocol and reconstruction (CASPAR) -reduction of image noise and patient dose[J]. IEEE transactions on medical imaging, 2010, 29(3):724-732.
[2] JIANG H Q, LI W X, LIU Y M, et al. Comparison study of filters for poisson noise removal[C]∥2011 International Workshop on Nonlinear Circuits, Communication and Signal Processing. Tianjin office of international cooperation: Tanjin Uniterity, 2011:171-174.
[3] FANG Y, ZHOU Y B, GE D W. De-noising based on wavelet analysis and bayesian estimation for low-dose X-ray CT [C]∥Proceedings 2009 9th International Conference on Electronic Measurement & Instruments. 北京：中国电子学会，2009:526-1051.
[4] MAKITALO M, FOI A. Optimal inversion of the anscombe transformation in low-count poisson image denoising [J]. IEEE transactions on image processing, 2011, 20(1):99-109.
[5] JIANG H Q, ZHANG Y Y, MA L, et al. A shearlet-based filter for low-dose mammography [J]. Breast imaging, 2014,8539:707-714.
[6] 张云逸，蒋慧琴，马岭. 基于剪切波的自适应泊松噪声的去除算法[C]∥第十七届全国图象图形学学术会议.珠海：中国图像图形学会，2014：395-401.
[7] ZHANG A, JIANG H Q, MA L, et al. A shearlet-based algorithm for quantum noise removal in low-dose CT images[J]. Proceedings of the SPIE, 2016, 9784: 978430.
[8] DAN Z P, CHEN X, GAN H T, et al. Locally adaptive shearlet denoising based on bayesian MAP estimate [C]∥IEEE Sixth International Conference on Image and Graphics. Hefei: IEEE Press, 2011:28-32.
[9] 胡海智,孙辉,邓承志,等.全变差正则化的Shearlet收缩去噪[J].中国图像图形学报, 2011, 16(2):168-173.
[10] LUISIER F, VONESCH C, BLU T, et al. Fast interscale wavelet denoising of poisson-corrupted images[J]. Signal process, 2010, 90(2): 415-427.
[11] MAKITALO M, FOI A. Optimal inversion of the anscombe transformation in low-count poisson image denoising[J]. IEEE transactions on image processing, 2011, 20(1):99-109.
[12] ZHANG B, FADILI J M, STARCK J L. Wavelets, ridgelets, and curvelets for poisson noise removal[J].IEEE transactions on image processing， 2008, 17(7): 1093-1108.