张　芳　桂志国,2*　张　权　董婵婵

1(中北大学信息与通信工程学院电子测试技术国家重点实验室　山西 太原 030051)2(中北大学信息与通信工程学院仪器科学与动态测试教育部重点实验室　山西 太原 030051)



基于小波和四阶各向异性扩散的MLEM低剂量CT重建算法

张芳1桂志国1,2*张权1董婵婵1

1(中北大学信息与通信工程学院电子测试技术国家重点实验室山西 太原 030051)2(中北大学信息与通信工程学院仪器科学与动态测试教育部重点实验室山西 太原 030051)

在低剂量计算机断层扫描CT(computedtomography)重建算法中，传统的最大似然期望最大MLEM(MaximumLikelihoodExpectationMaximization)算法随着迭代次数的增加会出现棋盘效应而不能有效地抑制噪声。针对上述问题提出一种基于小波收缩和四阶各向异性扩散相结合的MLEM低剂量CT重建算法。该算法结合小波收缩和各向异性扩散的优点，在每次迭代中，对MLEM重建算法处理后的图像进行离散平稳小波分解，在小波域的高频部分进行小波收缩，低频部分使用降噪效果优质的四阶各向异性扩散进行消噪，最后残留的脉冲噪声点通过中值滤波器进行处理，从而进一步优化图像。仿真实验结果表明，该算法可以有效地去除低剂量CT图像的噪声，且在保持图像边缘和细节信息方面有很好的表现，从而获得高抗噪性能的图像。

低剂量计算机断层扫描四阶各向异性扩散图像重建中值滤波

0　引　言

X射线计算机断层CT成像以其成像分辨率高的优点，在医学临床诊断中已获得广泛应用。然而，因其较高的辐射剂量会对人体造成一定的伤害，所以在扫描时应尽可能地降低射线剂量，即采用低剂量CT的成像技术。然而，低剂量CT重建会使图像质量发生明显的退化，为了解决上述问题,目前国内外学者尝试各种各样的处理方法来提高低剂量CT重建图像的质量，使其达到满意的效果。

针对低剂量CT重建图像的退化问题，目前主要方法可归于以下两种方式：1) 在重建后对得到的图像进行图像域降噪。2) 为先对图像的投影数据进行投影域降噪，然后再用降噪后的投影数据进行重建。近年来，针对前者方法学者们提出很多方法，且取得了令人满意的效果。如ChenYang[1]等通过使用一种新颖的非局部自适应加权非局部先验统计重建方法改善了低剂量CT图像的质量；Rust[2]等通过使用非线性高斯滤波器链对低剂量CT图像进行滤波，从而得到保护了边缘信息的降噪处理；LuiD[3]提出一种新颖的噪声补偿CT重建方法，从而提高了重建图像的信噪比。对于后者方法，即先对图像的投影数据进行投影域降噪，然后再用降噪后的投影数据进行重建，也取得了很大的成果。如由于传统的滤波反投影算法FBP(FilteredBackProjection)的算法简单，且速度较快，故该类算法一般是利用低剂量投影数据的特点设计滤波算法，滤除噪声后再利用FBP进行重建，从而获得良好的CT图像质量。GuiZhiguo[4]等使用模糊滤波器对低剂量CT的投影数据进行降噪处理。ZhangQuan[5]等对基于各向异性扩散加权先验的贝叶斯正弦图进行平滑。刘祎[6]等通过把局部模糊熵的自适应算法加入到投影数据中，从而得到高信噪比的重建图像。马建华[7]等通过对投影数据进行统计建模，然后采用贝叶斯最大后验估计的方法，将其非局部的先验信息加入到投影数据中，从而实现投影数据降噪的效果。WangJing[8]等利用惩罚加权最小二乘法方法，分别在图像域、投影域以及在投影数据的K-L(Karhunen-Loeve)域中进行了惩罚加权最小二乘法PWLS(PenalizedWeightedLeast-Squares)，取得了不错的效果。

虽然投影域降噪的重建方法取得了较大的发展，但由于在投影域进行降噪后，即使投影域中很小的点噪声映射到图像域，也会形成条形伪影。因此本文提出一种基于四阶各向异性扩散和小波收缩的的低剂量CT重建算法。该算法重建出来的图像不仅能够有效地保持图像边缘的细节信息，而且不存在各向异性扩散存在的明显阶梯效应。

1　本文重建算法

1.1MLEM重建方法

最大似然期望最大化MLEM算法由于在重建过程中既考虑了系统的物理模型，又考虑了观测数据的统计特性，因而被广泛的使用，其重建公式如下式：

(1)

1.2重建过程中的降噪方法

1.2.1四阶偏微分方程去噪算法

基于偏微分方程PDE(PartialDifferentialEquation)的各向异性扩散的去噪方法，可以满足整幅图像中需要去噪的不同强度的需求，因此是一种自适应的去噪技术。该方法的原理是在图像的平滑区域增加平滑的强度，而在图像的边缘位置则适当地减弱平滑强度从而保持图像的边缘区域，从而可以更好地在滤除噪声的同时抑制图像边缘的过平滑。由于二阶PDE降噪会出现“阶梯”效应，因此后人提出了四阶PDE模型[9]，其表达式为：

(2)

其中，fηη和fξξ为图像梯度方向和切线方向的二阶导数，表达式如下：

(3)

(4)

式(2)的离散化形式为：

(5)

四阶PDE既可以克服二阶PDE产生的“阶梯”效应，又可以根据梯度和切线方向进行不同程度的扩散，进而有效地保持边缘。

1.2.2小波变换

目前，图像降噪的方法有很多，小波变换作为一种新的多分辨分析方法，在图像去噪领域被广泛应用。小波变换能将信号的能量集中到少数小波系数上，而白噪声在任何正交基上的变换仍然是白噪声，并且有着相同的幅度。文献[10]分析小波系数在层间存在较强的持续性，小波域中信号系数随着尺度增加而增加，而噪声系数随着尺度的增加而减小。由于噪声和图像边缘主要集中在图像的高频部分，因此可以采用小波收缩法达到图像去噪的目的。

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

(6)

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

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

(7)

其中，wf(j,k)、ws(j,k)和wn(j,k)分别表示包含噪声信号、原始信号和噪声信号在第j层上的小波系数。由式(7)可以看出，观测信号的小波系数wf(j,k)，由原始信号小波系数ws(j,k)和噪声信号小波系数wn(j,k)两部分组成。对包含噪声信号的小波系数，若它的值大于指定的阈值，则此系数为信号分量，予以保留；若其值小于阈值，则认为此系数为噪声分量，滤除这样的系数就可达到降噪效果。

上述方法为小波阈值去噪，又称小波收缩[11]，它对小波系数进行统一处理，在最小均方误差下达到近似最优。和正交小波变换相比，平稳小波变换具有平移不变性而更适合用于图像的处理[12，13]。因此本文选用平稳小波变换对图像进行处理。其中2j为相应的尺度。

1.3本文提出的重建算法

由于小波变换具有良好的时频局部特性，小波变换之后，低频部分所含的噪声比较少。高频部分主要包含了图像的噪声和边缘，且噪声对应的小波系数幅值在噪声处比较小，在边缘处比较大。四阶偏微分对噪声的敏感性高，噪声越少，降噪效果越好。小波阈值收缩对处理此处高频噪声有很好的效果，故把四阶偏微分降噪用在低频处理，小波阈值用在高频处理是合适的。根据上述分析，本文的流程图1所示。

图1　本文算法流程图

本文的具体重建算法如下：

1)MLEM重建算法：

(8)

2) 小波收缩和四阶各向异性降噪处理：

在每次重建迭代中，首先对上步重建后的图像信号进行多尺度小波变换，生成相应的低频分量CAi，高频分量CHi,CVi,CDi，i=1,2,3,…,n为分解尺度。然后，对高频系数进行软阈值处理，去除噪声；对低频系数用基于四阶各向异性扩散算法进行图像降噪处理。所需公式如下：

(9)

(10)

其中，式(9)为软阈值收缩函数，ω表示含噪的小波系数，δ(ω)表示去噪后的小波系数，T为阈值，经过与其他几种常用的阈值做对比试验，证明本文选取的长度对数阈值效果最优且能达到令人满意的效果。故本文使用长度对数阈值，即对小波分解的高频小波系数按式(9)进行处理，对低频部分按式(10)进行四阶各向异性扩散处理，最后进行平稳小波逆变换重构，得到优化后的离散图像。

3) 进一步中值滤波处理：

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

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

(11)

所取窗口太大，不仅会增加算法的运算量，而且会使图像过平滑，所取窗口太小，会达不到降噪的效果，经过反复试验，这里取w为3×3的窗口。

4) 重复以上步骤，反复调试，选取视觉效果最优的结果作为最终的重建图像。本文提出算法的部分代码如下：

fork=1:150

F=(F·*(ggg*(p./(G*F))));

g=reshape(F,128,128);

%对MLEM重建后的图像进行平稳离散小波变换

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

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

sorh= ′s′;

%对高频成分进行软阈值处理

dswh=wthresh(swh,sorh,thr);

dswv=wthresh(swv,sorh,thr);

dswd=wthresh(swd,sorh,thr);

h=0.8;

%对低频成分进行四阶各向异性降噪处理

fork1=1:3

g=swa(:,:,k1);

k=4;

fori=2:M-1

forj=2:N-1

gx=(g(i+1,j)-g(i-1,j))/2;

gy=(g(i,j+1)-g(i,j-1))/2;

gxx=g(i+1,j)+g(i-1,j)-2*g(i,j);

gyy=g(i,j+1)+g(i,j-1)-2*g(i,j);

gxy=(g(i+1,j+1)+g(i-1,j-1)-g(i+1,j-1)-g(i-1,j+1))/4;

tidu=sqrt(gx*gx+gy*gy);

gc=(gx*gx*gxx+2*gy*gx*gxy+gy*gy*gyy)/(tidu*tidu+eps);

gs=(gy*gy*gxx-2*gy*gx*gxy+gx*gx*gyy)/(tidu*tidu+eps);

C=(k*k)/(k*k+(tidu*tidu));

l(i,j)=C*C*gc+C*gs;

end

end

fori=2:M-2

forj=2:N-2

t(i,j)=l(i+1,j)+l(i-1,j)+l(i,j+1)+l(i,j-1)-4*l(i,j);

g(i,j)=g(i,j)-dt*t(i,j);

end

end

swa(:,:,k1)=g;

end

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

2　实验结果与分析

2.1重建图像比较

本文算法及所有比较算法的实验仿真环境为：Windows7 旗舰版32位SP1(DirectX11)操作系统，英特尔Celeron(赛扬)E3300@2.50GHz双核处理器，2GB内存。编程工具使用MATLAB7.6.0(R2008a)。本文选取大小为128mm×128mm的Sheep-Logan头部剖面图模型作为实验对象，其灰度范围为[0,255]，如图2(a)所示。实验的无噪声投影数据通过在180个角度中均匀选取128个投影角度，每个投影角度下有128条射线的方法来获得，故本文选取大小为(128×128)×(128×128)的系统矩阵，实验中采用平行投影方式，在180个角度中均匀采取128个投影方向，每个方向分布128个探测器对。本文按照下式向理想投影数据加入期望和方差成非线性关系的高斯噪声来仿真低剂量CT投影数据：

σi2=kiexp(λi/T)

(12)

其中，i=1,2，…，N表示探测器信道，N表示信道总数，λi表示第i个探测器获得投影数据的均值，σi2表示第i个探测器获得投影数据的方差，ki表示第i个探测器的参数，T表示用于描述扫描系统校准过程的系统参数。对于给定的CT采集系统，ki与T是给定的。ki越大，T越小，低剂量CT重建图像的噪声越大；反之ki越小，T越大，低剂量CT重建图像的噪声越小。若噪声太大，会使图像失真，从而得不到满意的效果；若噪声太小，不能达到和实际数据的比拟，失去了研究的实际意义。为了更好地模拟实际低剂量CT的临床数据，经过大量的实验经验与对比，本文噪声的参数选取：ki=200，T=12 000。为了验证本文算法的有效性，将本文算法与传统的MLEM、BI-PLS(blockiteration-PenalizedLeastSquares)[15]、OS-PML-OSL(OrderedSubsets-PenalizedMaximumLikelihood-OneStepLate)和在每次重建迭代中使用基于四阶的各向异性扩散降噪的算法进行比较。各种算法的对比结果如图2所示。

图2　各种算法的对比结果

图2中，(a)为原始图像，(b)采用传统的MLEM重建算法，(c) 采用基于块迭代的BI-PLS重建算法，(d)采用基于有序子集的OS-PML-OSL重建算法。为了把本文算法和降噪效果很好的各向异性扩散算法进行比较，(e)在每次迭代中使用基于四阶的各向异性扩散降噪算法对重建图像进行处理，(f)为本文提出的算法。各种算法中涉及到的各种参数以及迭代次数，均为反复实验后得到的最优值，从结果图中可得，传统的MLEM重建图像的质量最差；BI-PLS和OS-PML-OSL算法的结果图可以达到降噪的效果，但同时也模糊了图像的边缘和细节；基于四阶的各向异性扩散降噪算法对重建图像的噪声进行了一定程度的抑制且可获得比较清晰的图像，但是图像中存在一些比较明显的块状阴影；本文算法在消除噪声的同时很好地保持了图像的边缘和细节信息，使图像达到比较优质的效果，从视觉上分析，重建效果和其他几种比较算法相比达到最优，初步说明了本文算法的有效性。

2.2重建精度比较

从2.1节分析可知，本文提出的算法对低剂量CT的重建图像有很好的降噪以及保持边缘和细节信息的能力，为了进一步说明本文算法的有效性，本文采用归一化均方距离、均方绝对误差、信噪比对重建图像质量进行定量描述。其定义如下：

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

(13)

2) 均方绝对误差MAE(MeanAbsoluteError):

(14)

3) 信噪比SNR(SignaltoNoiseRatio):

(15)

其中，J表示图像像素点的总数，Fi和fi分别表示重建图像与原始图像的第i个像素的灰度值，Mi和mi分别表示重建图像与原始图像的平均值。这些指标分别从不同的方面评价重建图像与原始图像的接近程度以及重建图像的质量，表1为本文算法与其他几种比较算法的客观评价结果。

表1　各种算法的客观评价

从表1可得出本文算法的信噪比最大，其他指标值均比其他的几种比较算法小。该结论说明本文算法的重建图像和原始图像最为接近。因此在定量评价方面，同样表明本算法在低剂量CT重建中的可行性。

图3给出了本文所用的Sheep-Logan模型的原始理想图像与在以上各种算法的重建图像的侧面轮廓线的比较图。从图中可以看出本文算法的重建图像与原始图像的吻合度最高，最接近于理想图像，具有最小的噪声波动，可以在保持图像边缘和细节信息的同时较好地解决低剂量CT重建图像的噪声问题。

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

3　结　语

本文针对传统的MLEM重建算法收敛速度慢，不能很好地控制噪声的问题，提出了一种基于小波收缩和四阶各向异性扩散的MLEM的低剂量CT重建算法。其中小波收缩可以在去除低剂量图像噪声的同时很好地保持图像的边缘和细节信息；基于四阶的各向异性扩散可以有效地判断图像的背景和边缘，然后进行不同程度的扩散；在低剂量CT重建中，小噪声可以由上述方法处理，而大噪声等脉冲噪声则可以利用中值滤波器去除。本文算法结合了以上各种算法的优势，既有效地保持了图像的细节和边缘信息，又很好地解决了低剂量图像的噪声问题。实验结果表明，该算法在主观效果和客观效果来看，均说明本文算法是可行的。在临床中该算法可以提高低剂量CT重建图像的质量，使其达到满意的效果，使医生在患者接受较少的辐射剂量下便可进行良好的诊断和治疗。本算法不仅可以应用于医学CT的医学成像领域，也可用于工业CT领域，同时在医学PET和工业无损检测中也有良好的应用前景。

[1]ChenY,GaoDZ,NieC,etal.Bayesianstatisticalreconstructionforlow-doseX-raycomputedtomographyusinganadaptive-weightinglocalnonprior[J].ComputerizedMedicalImagingandGraphics,2009,33(7): 495-500.

[2]RustGF,AurichV,ReiserM.Noisedosereductionandimageimprovementsinscreeningvirtualcolonoscopywithtubecurrentsof20mAswithnonlinearGaussianfilterchains[C]//MedicalImaging2002Conference,NewYork:IEEE,2002:186-197.

[3]LuiD,CameronA,ModhafarA,etal.Low-dosecomputedtomographyviaspatiallyadaptiveMonte-Carloreconstruction[J].ComputerizedMedicalImagingandGraphics,2013,37(7-8): 438-449.

[4]GuiZG,LiuY.Noisereductionforlow-doseX-raycomputedtomographywithfuzzyfilter[J].Optik-InternationalJournalforLightandElectronOptics,2012,123(13):1207-1211.

[5]ZhangQ,GuiZG,ChenY,etal.Bayesiansinogramsmoothingwithananisotropicdiffusionweightedpriorforlow-doseX-raycomputedtomography[J].Optik-InternationalJournalforLightandElectronOptics,2013,124(17):2811-2816.

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

[7] 马建华，黄静，陈阳，等.基于广义Gibbs先验的低剂量X-CT优质重建研究[J].计算机工程与应用,2008,44(16):4-6,93.

[8]WangJ,LiTF,LuHB,etal.PenalizedWeightedLeast-SquaresApproachtoSinogramNoiseReductionandImageReconstructionforLow-DoseX-RayComputedTomography[J].IEEETransactionsonmedicalimaging,2006,25(10):1272-1283.

[9]HajiaboliMR.AnAnisotropicFourth-OrderDiffusionFilterforImageNoiseRemoval[J].InternationalJournalofComputerVision,2011,92(2):177-191

[10]LiuJ,MoulinP.Information-theoreticanalysisofinterscaleandintrascaledependenciesbetweenimagewavelet-coefficients[J].IEEETransImageProcessing,2001,10(11):1647-1658.

[11]LeeK,VidakovicB.Semi-supervisedwaveletshrinkage[J].ComputationalStatistics&DataAnalysis,2012，56(6)：1681-1691.

[12]DengJM,YueHZ,ZhuoZZ,etal.AstationarywavelettransformbasedapproachtoregistrationofplanningCTandsetupconebeam-CTimagesinradiotherapy[J].JournalofMedicalSystems,2014,38(5)：1-9.

[13]JumahAA.Denoisingofanimageusingdiscretestationarywavelettransformandvariousthresholdingtechniques[J].JournalofSignalandInformationProcessing,2013,4(1)：33-41.

[14]LingJ,BovikAC.Smoothinglow-SNRmolecularimagesviaanisotropicmedian-diffusion[J].IEEETransactionsonMedicalImaging,2002,21(4):377-384.

[15]ByrneCL.Block-Iterativemethodsforimagereconstructionfromprojections[J].IEEETransactionsonimageprocessing,1996,5(5):792-794.

MLEMLOW-DOSECTRECONSTRUCTIONBASEDONWAVELETANDFOURTH-ORDERANISOTROPICDIFFUSION

ZhangFang1GuiZhiguo1,2*ZhangQuan1DongChanchan1

1(National Key Laboratory for Electronic Measurement Technology,School of Information and Communication Engineering,North University of China,Taiyuan 030051,Shanxi,China)2(Key Laboratory of Instrumentation Science and Dynamic Measurement,School of Information and Communication Engineering,North University of China,Taiyuan 030051,Shanxi,China)

Inlow-doseCT(computedtomography)reconstructionalgorithm,traditionalmaximumlikelihoodexpectationmaximisation(MLEM)algorithmwillappearchessboardeffectalongwiththeincreaseofthenumberofiterations,thuscannoteffectivelysuppressnoises.Forthisproblem,thispaperproposesalowdoseCTreconstructionalgorithm,whichisbasedonthecombinationofwaveletshrinkageandfourth-orderanisotropicdiffusion.Itcombinestheadvantagesofwaveletshrinkageandanisotropicdiffusion,ineachiteration,itconductsthediscretestationarywaveletdecompositionontheimageprocessedwithMLEMreconstructionalgorithm,inhighfrequencypartofthewaveletdomainitshrinksthewavelet,inlowfrequencypartitusesfourth-orderanisotropicdiffusion,whichhashighqualityeffectindenoising,toeliminatenoises,itprocessesthefinalresidualpulsenoisepointswiththemedianfilter,sothatfurtheroptimisestheimage.Simulationexperimentresultsshowthattheproposedalgorithmcaneffectivelyremovethenoiseinlow-doseCTimage,andhasgoodperformancesinkeepingboththeimageedgesanddetailedinformation,therebygainstheimagewithhighanti-noiseperformance.

Low-dosecomputedtomographyFourth-orderAnisotropicdiffusionImagereconstructionMedianfilter

2014-09-03。国家自然科学基金项目(61071192，61271357，61171178)；山西省国际合作项目(2013081035)；山西省研究生优秀创新项目(2009011020-2，20123098)；中北大学第十届研究生科技基金项目(20131035)；山西省高等学校优秀青年学术带头人支持计划资助项目；中北大学2013年校科学基金计划。张芳，硕士生，主研领域：基于低剂量CT的图像重建。桂志国，教授。张权，副教授。董婵婵，硕士生。

TP391

ADOI:10.3969/j.issn.1000-386x.2016.03.043