董婵婵　桂志国　张　权　郝慧艳　张　芳　刘　祎　孙未雅

(中北大学信息与通讯工程学院　山西 太原 030051)



基于ROAD和小波收缩的MLEM低剂量CT重建算法

董婵婵桂志国张权郝慧艳张芳刘祎孙未雅

(中北大学信息与通讯工程学院山西 太原 030051)

摘要针对低剂量CT (Computed Tomography)重建图像质量退化的问题，提出一种基于小波收缩和绝对差值排序各项异性扩散的MLEM(Maximum Likelihood Expectation Maximization)低剂量CT重建算法。算法在每次迭代中首先采用MLEM算法对低剂量CT投影数据进行重建。由于各项异性扩散对噪声敏感，所以算法先对重建后的图像进行小波变换，再在更稳定的低频小波域进行基于绝对差值排序的各项异性扩散处理，对小波高频系数进行软阈值降噪处理。然后将降噪处理后的系数进行小波反变换，得到降噪后的图像。最后使用中值滤波对图像进行处理，从而消除脉冲噪声点。实验结果表明，与其他几种常用重建算法相比，该算法重建的图像信噪比更高，归一化均方误差更小，处理后的图像更清晰，即可以在抑制噪声的同时，较好地保持图像的边缘和细节信息。

关键词低剂量CT小波变换图像重建绝对差值排序各向异性扩散

MLEM LOW-DOSE CT RECONSTRUCTION ALGORITHM BASED ON ROAD AND WAVELET SHRINKAGE

Dong ChanchanGui ZhiguoZhang QuanHao HuiyanZhang FangLiu YiSun Weiya

(School of Information and Communication Engineering, North University of China, Taiyuan 030051,Shanxi,China)

AbstractConcerning the problem of quality degradation of low-dose CT reconstruction images, we presented an MLEM low-dose CT reconstruction method which is based on wavelet shrinkage and rank-ordered absolute differences anisotropic diffusion. In each time of iteration, the algorithm first uses MLEM to reconstruct the low-dose projection data. Since the anisotropic diffusion is sensitive to noises, so the algorithm performs wavelet transform on the reconstructed image prior to conducting anisotropic diffusion processing based on rank-ordered absolute differences in more stable low-frequency wavelet domain and then carries out the soft threshold denoising processing on high-frequency coefficient of wavelet. After that the algorithm performs inverse discrete wavelet transform (IDWT) on the coefficients with denoising treatment and obtains the denoised images. Finally it uses median filter to process the image so as to eliminate the impulse noise points. Experimental results showed that compared with some other common reconstruction algorithms, the image reconstructed by the proposed one had higher signal-to-noise ratio (SNR) and smaller normalised mean square error, the processed image with the proposed algorithm were much clear, i.e., it could well preserve edges and details information of image while suppressing the noise.

KeywordsLow-dose computed tomography (CT)Wavelet transformImage reconstructionRank-ordered absolute differences (ROAD)Anisotropic diffusion

0引言

计算机断层成像术CT技术在疾病预防、临床诊断中发挥了巨大的作用[1]。然而，因其较高的辐射剂量会对人体造成一定的伤害，所以在扫描时应尽可能地降低射线剂量，即采用低剂量CT成像技术[2]。低剂量CT可减少对受检者的放射危害，然而当剂量降低时，由于穿过人体抵达探测器的光子数目过少，投影数据会被随机噪声严重污染，导致重建得到的图像质量退化明显，重建图像中可见明显的条状伪迹。因而，有效地滤除低剂量CT图像的噪声，有着重要的实际应用价值。

到目前为止，人们已经提出了很多滤波策略来解决噪声问题，主要包括在投影域进行统计滤波并用滤波反投影算法重建图像和在CT重建图像域直接去噪的算法。目前在投影域进行滤波的研究工作主要包括：Kachelriess M等人[3]提出了一种投影域的自适应卷积滤波方法；毕一鸣等人利用噪声的统计特性[4]，提出了一种基于标准剂量CT图像非局部权值先验的低剂量图像恢复方法； Gui Zhiguo[5]等根据低剂量CT噪声的特性，提出一种模糊中值滤波器，用此对低剂量CT的投影数据进行降噪处理，从而提高图像的抗噪声性能；刘祎[6]等通过使用局部模糊熵判断细节和平滑区域，然后利用基于统计信息的各向异性滤波作用在投影数据中，从而得到高信噪比的恢复图像。由于在投影域降噪后，即使遗留下很小的点噪声，映射到图像域也会形成条形伪影。因此，在重建后的图像上进行去噪的算法得到了广泛关注。Chen Yang[7]等通过使用非局部自适应加权先验统计方法进行图像重建，改善了低剂量CT图像的质量。王丽艳[8]等通过采用待重建图像稀疏性的先验信息作为正则项，泊松噪声的负对数似然函数作为保真项设计优化目标函数，从而达到去噪且保持细节的目的。Rust 等利用非线性高斯滤波器链对重建图像进行平滑处理[9]，得到不错的降噪效果且很好地保持了边缘和细节。Lui D[10]提出一种新颖的噪声补偿CT重建方法从而提高了重建图像的信噪比。

目前图像降噪的方法有很多，如：数学形态学法、小波技术[11，12]、非局部均值滤波法[13]以及基于偏微分方程的图像降噪算法法[14，15]。近年来，基于小波变换和基于偏微分方程的图像处理方法受到人们的广泛关注。通过以上研究分析，本文提出了一种基于小波和绝对差值排序ROAD(rank-ordered absolute differences)各项异性扩散的最大似然期望最大化MLEM的低剂量CT重建算法。该算法先对投影数据进行重建，在将重建得到的图像进行小波分解，对高频系数进行阈值降噪处理，对小波低频系数进行基于ROAD的各项异性扩散。实验结果表明，该算法重建出来的图像可以在对图像进行降噪的同时，较好地保持图像的边缘和细节信息。

1噪声模型

关于投影数据噪声模型的数据特点已经进行了广泛的研究，文献[16]认为探测到的光子数服从背景噪声为高斯噪声的泊松分布。Lu等人通过对多个体模进行反复扫描获得的投影数据进行分析得出[17，18]：经系统校准及对数变换之后的低剂量CT投影数据近似服从非平稳高斯分布，且其均值和方差之间具有非线性关系。其统计模型为：

σi2=kiexp(λi/T)

(1)

其中，i=1,2,…,N表示探测器信道，N表示信道总数，λi表示第i个探测器获得投影数据的均值，σi2表示第i个探测器获得投影数据的方差，ki表示第i个探测器的参数，T表示用于描述扫描系统校准过程的系统参数。对于给定的CT采集系统，ki与T是给定的。

2基于绝对差值排序和小波收缩的MLEM低剂量CT重建算法

2.1MLEM重建方法

MLEM算法由于在重建过程中考虑了观测数据的统计特性，使重建出的图像优于FBP重建出的图像，重建公式如下：

(2)

2.2本文算法的理论知识

2.2.1基于绝对差值排序的各项异性扩散

P-M方程的表达式为：

(3)

其中ct为扩散系数函数，其表达式为：

(4)

由于有些噪声的梯度值与边缘的梯度值相差无几，甚至更大，因此只依靠梯度信息来检测边缘是不够的。文献[15]指出：即使在噪声与边缘的梯度值相似时，噪声与边缘处的灰度方差却相差很大，细节处的灰度方差明显比噪声处的灰度方差大。Chao和Tsai在P-M模型的基础上增加了灰度方差的信息，可以在去除噪声的同时保留更多的细节。

通过在扩散过程中包含梯度和灰度方差的局部特征，式(4)中的扩散系数函数改为：

(5)

为了表示出某个特定的像素值与它周围的邻域像素值强度的差异，Roman Garnett等人[19]提出了一个局部图像统计—绝对差值排序检测法ROAD。设中心像素点的坐标为x=(x1,x2),则：

Ωx(N):={x+(i,j):-N≤i,j≤N}

(6)

(7)

由于图像的内部区域及边缘都具有连续性，所以在它周围的八邻域中，至少有四个邻域值与中心像素的强度相似，说明其ROAD值较小。而噪声会使中心像素与大多数或所有的邻域像素在强度上相差很大，说明ROAD值较大。因此，通过ROAD值的大小可以区分出边缘与噪声。此时，扩散系数函数的表达式为：

(8)

基于绝对差值排序的各项异性扩散方程的离散化形式为：

ft+1(x,y)=ft(x,y)+

(9)

2.2.2小波降噪

图像降噪的目的是为了改善图像质量，减少噪声，同时尽可能保持图像的特征不被破坏。由于噪声和图像的特征都是高频成分，因此在降噪的同时会损失图像的细节特征。近年来，在图像去噪领域，小波理论受到了许多学者的重视，将小波应用到降噪领域，获得了非常好的效果。其中，小波阈值收缩法是目前研究最为广泛的小波去噪算法。由于小波系数在层间存在较强的持续性，小波域中信号系数随着尺度的增加而增加[20]，而噪声系数随着尺度的增加而减小。因此将信号进行小波分解后，高频小波系数中绝对幅值比较大的小波系数主要由真是图像信号提供，予以保留；绝对幅值比较小的小波系数主要由噪声提供，滤除这样的系数就可达到降噪效果。

f(k)=s(k)+n(k)k=0,1,2,…，N-1

(10)

对f(k)做离散小波变换，可得：

wf(j,k)=ws(j,k)+wn(j,k)

(11)

其中，wf(j,k)、ws(j,k)和wn(j,k)分别表示包含噪声信号、原始信号和噪声信号在第j层上的小波系数。

2.3算法描述

一方面，图像经过小波变换后，边缘和噪声是高频分量，图像的主要信息都集中在低频部分。由于各项异性扩散对噪声较敏感，而低频分量中的噪声较少，因此在低频部分进行各项异性扩散可以减小噪声对其的影响。另一方面，小波收缩去噪方法是对小波系数进行统一处理，同时考虑了小波系数之间的相关性和局部特征。非线性扩散方法与像素点的具体位置有关[21]，它是一种从局部向整体逐渐扩散的过程。但是这种方法需要多次迭代，计算量较大，相比之下，小波收缩的速度较快。由于小波收缩方法收敛速度快，所以在高频分量部分采用小波收缩软阈值去噪；从计算复杂性考虑，由于图像经过一次小波变换后，低频部分大小减为原始图像的1/4。所以，图像经过多尺度分解后，低频部分大小远小于原图像大小，因此只对低频部分采用非线性滤波，比在整幅图中采用非线性滤波所用的时间大大减少。因此，本文算法的降噪部分采用混合图像去噪算法。

算法先对投影数据进行重建，然后将重建后的图像进行小波分解，使得图像的主要信息大都集中在低频，图像的边缘和噪声大都集中在高频。再对信号的低频小波系数进行各项异性扩散处理，减少噪声对各项异性扩散的影响；对小波域中信号的高频系数进行软阈值处理。将降噪处理后的小波系数进行小波反变换，最后再对降噪后的图像进行中值滤波处理。由于噪声和边缘都是高频信息，但在小波域中信号小波系数随着尺度的增加而变大，而噪声系数随着尺度的增加而减小，所以通过软阈值处理可以有效去除噪声。由于图像的信息主要集中在低频，低频中包含的噪声较少，因此本文将差分曲率算子引入到各项异性扩散方程中，对低频信息进行小幅度的平滑处理，去除低频中的噪声。根据上述分析，本文算法的流程如图1所示。

本文根据图1中的流程可知该算法的具体步骤如下：

(2) 在每次重建迭代中，首先对上步重建后的图像信号进行多尺度小波变换，生成相应的低频分量CAi，高频分量CHi,CVi,CDi，i=1,2,…,n为分解尺度。

(3) 然后，对高频系数进行软阈值处理，去除噪声；对低频系数用基于绝对差值排序检测各项异性扩散算法进行图像降噪处理。

阈值处理过程分为软阈值和硬阈值，硬阈值方法可以很好地保留图像边缘等局部特征，但图像会出现振铃、伪吉布斯效应等视觉失真，而用软阈值方法处理结果则相对平滑得多。因此，本文采用软阈值方法。

(4) 在小波域进行图像降噪处理后，进行小波反变换得到去噪后的图像。

(5) 对各向异性扩散处理后的脉冲噪声进行中值滤波处理：

fi,jn + 1=Median(fi,jn + 1,w)，

其中，w是中值算子的窗口。

图1　本文算法流程图

由文献[22]可知，低剂量CT重建图像的噪声还表现为一些脉冲噪声，各向异性扩散降噪技术对重建的图像进行降噪后，仅可以平滑图像的小梯度区域，而相对于周围区域的大梯度区域则保持不变，这些大梯度可能是边缘，也可能是图像的峰值噪声。而中值滤波器只会对大噪声峰值产生的大梯度起作用，边缘产生的大梯度将不会受到影响。因此在低剂量CT重建时，低噪声可以由各向异性扩散平滑，而大噪声等脉冲噪声则由中值滤波器所消除。

即采用各项异性扩散降噪算法对重建图像进行处理后，小噪声引起的小梯度区域被平滑，而大噪声峰值如脉冲噪声由该中值滤波进行滤除。由于中值滤波是将窗口内像素的灰度值进行排序，用排序所得的中间值代替窗口中心位置像素的灰度值，若窗口选择太大会使计算量过大，运行时间过长，不利于图像处理的实时性，因此这里w取3×3的窗口。

(6) 重复步骤(1)-(5)直至得到最终的重建图像。

先对投影数据进行MLEM重建，然后进行图像域降噪处理。降噪部分即先进行基于小波收缩和绝对差值排序各项异性扩散，然后进行中值滤波。

F=(F.*(ggg*(p2./(G*F))));

%MLEM重建方法

g=reshape(F,128,128); x=double(g); xx=x;

[swa,swh,swv,swd] = swt2(xx,3,′db1′);

%采用db1进行多尺度小波分解

[thr] = ddencmp(′den′,′wv′,xx);

%对高频部分进行小波软阈值处理

sorh = ′s′; dswh = wthresh(swh,sorh,thr);

dswv = wthresh(swv,sorh,thr);

dswd = wthresh(swd,sorh,thr);

h=0.8; k=swa(:,:,3);

%对低频部分进行基于ROAD的各项异性扩散降噪处理

scalval=0.085; kval=1; x=5;

he=fspecial(′gaussian′, 3, 9);

imgMean1=imfilter(k,he);

imgMean=imfilter(k,he);

for ii=2:N-1

for jj=2:M-1

for iii=ii-1:ii+1

for jjj=jj-1:jj+1

sigma(ii,jj)=(k(iii,jjj)-imgMean(ii,jj))*(k(iii,jjj)-imgMean(ii,jj));

end

end

sigma(ii,jj)=sigma(ii,jj)/9;

end

end

sigma=1+(sigma-min(sigma(:)))/(max(sigma(:))-min(sigma(:)))*254;

for ii=2:N-1

for jj=2:M-1

B=zeros(3,3);

a1=abs(k(ii-1,jj)-k(ii,jj));

a2=abs(k(ii,jj-1)-k(ii,jj));

a3=abs(k(ii+1,jj)-k(ii,jj));

a4=abs(k(ii,jj+1)-k(ii,jj));

a5=abs(k(ii+1,jj-1)-k(ii,jj));

a6=abs(k(ii+1,jj+1)-k(ii,jj));

a7=abs(k(ii-1,jj-1)-k(ii,jj));

a8=abs(k(ii-1,jj+1)-k(ii,jj));

A=[a1,a2,a3,a4,a5,a6,a7,a8];

B=sort(A);

ROAD=B(1)+B(2)+B(3)+B(4); kkk=0;

tidu=k(ii-1,jj)-k(ii,jj); tidu1=tidu; kk1=tidu*(x*ROAD*kval*kval)/((x*ROAD*kval*kval)+sigma(ii,jj)*sigma(ii,jj)*tidu1*tidu1+eps);

tidu=k(ii,jj+1)-k(ii,jj); tidu1=tidu; kk2=tidu*(x*ROAD*kval*kval)/((x*ROAD*kval*kval)+sigma(ii,jj)*sigma(ii,jj)*tidu1*tidu1+eps);

tidu=k(ii+1,jj)-k(ii,jj); tidu1=tidu; kk3=tidu*(x*ROAD*kval*kval)/((x*ROAD*kval*kval)+sigma(ii,jj)*sigma(ii,jj)*tidu1*tidu1+eps);

tidu=k(ii,jj-1)-k(ii,jj); tidu1=tidu; kk4=tidu*(x*ROAD*kval*kval)/((x*ROAD*kval*kval)+sigma(ii,jj)*sigma(ii,jj)*tidu1*tidu1+eps);

kkk=kk1+kk2+kk3+kk4;

k(ii,jj)=k(ii,jj)+kkk*scalval;

end

end

swa(:,:,3)=k;

xd=iswt2(swa,dswh,dswv,dswd,′db1′);

%小波反变换

g=reshape(xd,128,128);

g=medfilt2(g);

%中值滤波

F=reshape(g,16384,1);

3实验结果与分析

为了验证算法的有效性，本文采用Sheep-Logan体模作为实验模型进行低剂量CT图像重建的仿真，与BIPLS(Block Iterative-Penalized Least Squares)，OS-PML-OSL(Ordered Subsets-Penalized Maximum Likelihood-One Step Late)进行了比较，且为了进一步说明本文中改进各向异性扩散算法在MLEM的低剂量CT重建算法的有效性，将基于传统PM、基于方差应用到MLEM算法中，并和本文算法进行了比较。采用大小为128 mm×128 mm，灰度为0~255的Sheep-Logan头部体模。实验通过在360°范围内均匀选取984个投影角度，每个投影角度下有888条射线的方法来获得仿真的无噪声投影数据。在此基础上，按照式(1)的方法向理想投影数据中加入高斯噪声来模拟低剂量CT的投影数据，其中ki=200，T=12 000。实验仿真环境为：32位Microsoft Windows XP Professional 2002计算机操作系统，E5300@2.60 GHz英特尔奔腾双核处理器 ，2 GB内存。软件环境为：MATLAB R2008a。图2为不同算法重建的图像。

图2　各种算法的对比结果

图2(a)为原始图像。(b) 采用OS-PLS重建算法。(c) 采用OS-PML-OSL重建算法。(d) 使用MLEM算法进行重建，然后采用基于传统PM的各向异性扩散对图像域进行降噪。(e) 使用MLEM算法进行重建，然后采用基于方差的各向异性扩散对图像域进行降噪。(f) 采用本文提出的算法。

由图2(b)可以明显看出，MLEM重建算法的重建结果图中含有较多噪声，重建效果不好。(c)和(d)与(b)相比可以看出：算法BIPLS和OS-PML-OSL比(b) MLEM算法重建后的图像包含更少的噪声，重建效果较好。而(e)比(b)、(c)和(d)中包含更少的噪声，重建的图像更清晰。将(f)的结果与(b)、(c) 、(d) 和(e)相比可以明显看出：与其他四种算法的重建结果相比，本文提出算法的重建图像更清晰，包含的噪声更少。综上所述，本文算法可以有效地解决低剂量重建图像的噪声问题，且可以在光滑去噪的同时较好地保持图像的纹理和边缘信息。

图3(a)为一个大小为128 mm×128 mm的肩部仿真模型，(b)为MLEM重建后的结果，(c)为BIPLS算法处理后的结果，(d)为OS-PML-OSL处理后的图像，(e)为基于方差的MLEM处理后的结果，(f)为本文算法处理后的结果。

由图3可以清晰看出，图像(f)较清晰，图像质量较高，即本文处理后的图像明显优于其他算法处理后的图像。因此，本文算法可以有效地解决低剂量重建图像的噪声问题，且可以在光滑去噪的同时较好地保持图像的纹理和边缘信息。总之，从视觉上来看，重建效果优于其他几种比较算法。

图3　各种算法的对比结果

从上述分析可得，本文提出的算法在去噪能力与保持边缘和细节能力方面都明显优于其他重建方法。为了进一步证明本文算法的有效性，本文采用归一化均方距离、归一化均方误差、均方绝对误差、信噪比对其进行定量描述。其定义分别如下：

(1) 归一化均方距离NMSD:

(12)

(2) 归一化均方误差NMSE:

(13)

(3) 均方绝对误差MAE:

(14)

(4) 信噪比SNR:

(15)

其中，M和N分别为图像的行数和列数，Fi和fi分别表示重建图像与原始图像的第i个像素的灰度值，Mi和mi分别表示重建图像与原始图像的均值，本文提出的新算法以及各种比较算法关于Sheep-Logan体模仿真的客观评价结果如表1所示。

表1　各种算法的客观评价

由表1可以明显看出，本文提出算法的归一化均方误差NMSE、归一化均方距离NMSD、均方绝对误差MAE均比其他的几种比较算法小，信噪比SNR高达10.344450，比其他几种算法高出很多。因此由图2和表1可知，无论在视觉方面还是在定量评价方面，均表明本算法在CT重建中是可行且有效的。

图4给出了Sheep-Logan体模的以上重建图像与原始理想图像的侧面轮廓线的比较图。从图中可以看出本文算法对低剂量进行CT重建更为接近理想图像，与原始图像的吻合度最高，即和其他几种比较算法相比较，本文方法重建出的图像噪声波动最小，可以较好地去除低剂量CT图像的噪声以及较好地保持图像的细节和纹理信息。

图4　各种算法第65行侧面轮廓线的对比结果

4结语

本文提出了一种基于绝对差值排序各项异性扩散和小波收缩的MLEM的低剂量CT重建算法。各项异性扩散虽然可以有效去除图像的噪声，但它对噪声敏感。所以，算法在MLEM重建之后，先进行小波变换，然后在更稳定的小波域而不是噪声图像域进行基于绝对差分排序的各项异性扩散处理，减少噪声对各项异性扩散的影响。实验验证表明：与其他常用算法相比，该算法重建出的图像更干净，含有更少的噪声，信噪比更高，保持更多的图像细节，重建效果更好。该算法可以应用于医学CT领域中，能够在降低辐射剂量的同时，有效地重建出符合要求的图像。因此，该算法在医学成像领域具有较好的应用前景。

参考文献

[1] 张元科, 张军英, 卢虹冰.低剂量CT图像模型参数估计及统计去噪研究[J].西安电子科技大学学报,2011,38(3)：99-106.

[2] 张权，罗立民，桂志国.基于改进非局部先验的Bayesian低剂量CT投影平滑算法[J].东南大学学报，2014，44(3)：499-503.

[3] Kachelriess M, Watzke O, Kalender W A. Generalized multi-dimensional adaptive filtering for conventional and spiral single-slice, multi-slice, and cone-beam CT[J]. Medical Physics,2001,28(4):475-490.

[4] 毕一鸣, 马建华, 陈武凡.基于标准剂量CT图像非局部权值先验的低剂量图像恢复[J].电子学报,2010,38(5):1146-1151.

[5] Gui Z G, Liu Y. Noise reduction for low-dose X-ray computed tomography with fuzzy filter[J].Optik-International Journal for Light and Electron Optics, 2012,123(13):1207-1211.

[6] 刘祎, 张权, 桂志国.基于模糊熵的低剂量CT投影降噪算法研究[J].电子与信息学报, 2013,35(6):1421-1427.

[7] Chen Y, Gao D Z, Nie C, et al. Bayesian statistical reconstruction for low-dose X-ray computed tomography using an adaptive-weighting local nonprior[J]. Computerized Medical Imaging and Graphics, 2009,33(7):495-500.

[8] 王丽艳, 韦志辉. 低剂量CT 的线性Bregman 迭代重建算法[J].电子与信息学报,2013,35(10):2418-2424.

[9] Rust G F, Aurich V, Reiser M. Noise dose reduction and image improvements in screening virtual colonoscopy with tube currents of 20 mAs with nonlinear Gaussian filter chains[C]//Medical Imaging 2002 Conference, New York: IEEE,2002:186-197.

[10] Lui D, Cameron A, Modhafar A, et al. Low-dose computed tomography via spatially adaptive Monte-Carlo reconstruction[J].Computerized Medical Imaging and Graphics, 2013,37(7-8):438-449.

[11] 王凯, 刘甲甲, 袁建英,等.基于小波子带系数能量的优化权值降噪[J].计算机应用,2013,33(8):2341-2345.

[12] Zhong J, Sun H. Wavelet-based multiscale anisotropic diffusion with adaptive statistical analysis for image restoration[J].IEEE Transactions on Circuits and Systems,2008,55(9):2716-2725.

[13] 钟莹, 杨学志,唐益明,等.采用结构自适应块匹配的非局部均值去噪算法[J].电子与信息学报,2013,35(12):2908-2915.

[14] Perona P, Malik J. Scale space and edge detection using anisotropic diffusion[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,1990,12(7):629-639.

[15] Chao S M, Tsai D M. An improved anisotropic diffusion model for detail and edge preserving smoothing[J].Pattern Recognition Letters,2010,31(13):2012-2023.

[16] Elbakri I A, Fessler J A. Statistical image reconstruction for polyenergetic X-ray computed tomography[J].IEEE Transactions on Medical Imaging,2002,21(2):89-99.

[17] Lu Hongbing, Li Xiang, Li Lihong, et al. Adaptive noise reduction toward low-dose computed tomography[C].SPIE Proceedings Medical Imaging, 2003: Physics of Medical Imaging,2003,5030:759-766.

[18] Li Tian-fang, Li Xiang, Wang Jing, et al. Nonlinear sonogram smoothing for low-dose X-ray CT[J].IEEE Transactions on Nuclear Science,2004,51(5):2505-2513.

[19] Garnett R, Huegerich T, Chui C, et al. A universal noise removal algorithm with an impulse detector[J].IEEE Transactions on Image Processing,2005 ,14(11):1747-1754.

[20] Liu J, Moulin P. Information-theoretic analysis of interscale and intrascale dependencies between image wavelet-coefficients[J].IEEE Trans Image Processing,2001,10(11):1647-1658.

[21] 王知强. 基于小波收缩与非线性扩散的去噪算法[J].计算机工程, 2011,37(7):249-252.

[22] Ling J, Bovik A C. Smoothing Low-SNR Molecular Images Via Anisotropic Median-Diffusion[J].IEEE Transactions on Medical Imaging,2002,21(4):377-384.

中图分类号TP391

文献标识码A

DOI:10.3969/j.issn.1000-386x.2016.02.038

收稿日期：2014-08-07。国家自然科学基金项目(61071192,61271357,61171178)；山西省国际合作项目(2013081035)；山西省研究生优秀创新项目(2009011020-2)；山西省研究生优秀创新项目(20123098)；中北大学第十届研究生科技基金项目(20131035)；山西省高等学校优秀青年学术带头人支持计划项目；中北大学2013年校科学基金计划。董婵婵，硕士生，主研领域：图像降噪和图像重建。桂志国，教授。张权，副教授。郝慧艳，副教授。张芳，硕士生。刘袆，博士。孙未雅，硕士生。