基于改进各项异性扩散和小波变换的MLEM低剂量CT重建算法

董婵婵，张 权，郝慧艳，刘 祎

（中北大学 电子测试技术国家重点实验室，山西 太原030051）

0 引 言

医学计算机断层成像术（Computerized Tomography，CT）技术是通过测量物体不同角度下的射线投影，从而获得物体横截面信息的成像技术，被广泛应用于医疗诊断、无损检测等领域［1］.X射线的放射剂量越高，重建图像的质量越好.但在进行CT检查时，高剂量的辐射会对人体的健康造成危害；辐射剂量降低时，获得的投影数据会伴随着比较高的电子噪声，重建出的图像质量往往不能满足图像处理和医疗诊断的需求［2］.因此，在低剂量下重建出高分辨率和低噪声的CT图像的研究越来越受到关注.

在低剂量CT图像重建过程中，由于探测器接受到的光子数较少，投影数据受噪声污染严重，导致重建后图像的质量较差［3］.研究至今，在对低剂量CT重建图像进行降噪时，主要方法有两类：一是先对图像的投影数据进行投影域降噪，然后再用降噪后的投影数据进行重建；二是先重建，再对重建后得到的图像进行图像域降噪处理.到目前为止，针对投影域数据进行降噪的研究取得了很多成果.如Wang Jing等利用惩罚加权最小二乘法方法［4］，分别在图像域、投影域以及在投影数据的K-L（Karhunen-Loeve）域中进行了惩罚加权最小二乘法（Penalized Weighted Least-Squares，PWLS），取得了不错的效果；2008年，Lu Hongbing等人［5］在小波域中进行PWLS滤波算法，通过小波变换来对投影数据进行处理，大大提高了重建图像的质量；李凯旋［6］等提出一种基于投影数据恢复导引的双边滤波权值优化方法，在对噪声去除的同时有效的保护了图像的细节；Zhang Quan［7］等提出一种新颖的基于各项异性扩散加权先验，通过对贝叶斯统计的正弦图进行平滑，从而保持了图像的边缘.但对投影域降噪处理的方法计算复杂且计算量大，同时重建算法需要对系统的噪声、探测器与X线源间校准等过程进行准确重建，增加了优化难度，而优化重建算法难度较高不利于在高速CT系统中应用［8］.而在图像域进行降噪处理的方法比较方便直接，可以根据图像不同的特性提出不同的降噪方法.因此，近年来在图像域进行降噪处理成为了研究热点，并取得了较好的效果.如Rust等利用非线性高斯滤波器链对重建图像进行平滑处理［9］；Chen Yang［10］等通过使用一种新颖的非局部自适应加权非局部先验统计重建方法，改善了低剂量CT图像的质量.

针对以上研究分析，本文提出了一种基于小波收缩和差分曲率各项异性扩散的最大似然期望最大化（Maximum Likelihood Expectation Maximization，MLEM）的低剂量CT重建算法，先对投影数据进行重建，再将重建得到的图像进行小波分解，对高频系数进行阈值降噪处理，对小波低频系数进行基于差分曲率的各项异性扩散.

1 噪声模型

低剂量CT图像的噪声表现在图像上是一些噪声脉冲和条状伪影，从而使图像的信噪比下降严重，以及图像中的一些细节被污染而无法显现.低剂量投影数据的噪声特点是低剂量CT图像重建解决以上问题的关键.长期以来，投影数据的噪声都被认为符合泊松分布，这是基于光子数目符合泊松分布而确定的.而Li［11］等人经过多次研究分析认为低剂量的投影数据通过系统校正以及对数变换之后，其均值与方差之间是非线性关系，近似服从非平稳高斯分布.其满足如下公式

式中：i=1，2，…，M，表示的是探测器的信道，M为信道的总数；λi和σ2i分别表示探测器信道处获得的投影数据的均值和方差；fi表示信道探测器的参数；T是系统参数，它是一个用于描述系统校准过程的尺度系数.不同的CT采集系统，fi与T这两个参数也不相同，但是对于一个给定的CT采集系统，fi与T是已知的.

2 基于小波变换和各项异性扩散的MLEM低剂量CT重建算法

2.1 MLEM重建方法

MLEM算法由于在重建过程中既考虑了系统的物理模型，又考虑了观测数据的统计特性，重建出的图像要优于FBP重建出的图像，其重建公式为［12］

2.2 重建过程中的降噪方法

2.2.1 各项异性扩散算法

传统P-M模型的表达式为［13］

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

式中：K为常数.传统P-M模型仅利用梯度来控制模型的扩散速度，虽然可以有效地去除噪声，但图像中有些平坦区域的梯度与细节处的梯度相差无几，因此只用梯度来检测边缘是不够的.2.2.2 基于差分曲率的各项异性扩散

文献［14］提出的差分曲率可以较好地区分细节和平坦区域、独立噪声，其表达式为

式中：fηη和fξξ分别表示图像f在梯度方向和水平方向的二阶导数，其表达式为

由式（5）可知：在边缘处，|fηη|较大，|fξξ|较小，其差分曲率D值较大；在平坦区域，|fηη|和|fξξ|都较小，所以D值较小；在独立噪声点处，|fηη|和|fξξ|都较大，且几乎相等，所以D较小.因此根据差分曲率D值大小，可以很好地区分边缘和平坦区域、独立噪声点.因此，将差分曲率引入后的扩散系数函数为

其中，dt，N为归一化差分曲率

此时，各项异性扩散方程为

2.2.3 小波阈值降噪

小波变换是用小波系数的形式来表示信号，描述不同尺度的信号变化情况.由于小波变换具有时频局部化和多分辨率的特性，因此图像去噪在小波域中的研究成为热点.小波变换能将信号的能量集中到少数小波系数上［15-16］，而白噪声在任何正交基上的变换仍然是白噪声，并且有着相同的幅度.相对而言，信号的小波系数值必然大于那些能量分散且幅值较小的噪声的小波系数值.文献［17］中分析了小波系数在层间存在较强的持续性，小波域中信号系数随着尺度增加而增加，而噪声系数随着尺度的增加而减小.由于噪声主要集中在图像的高频部分，因此选择一个合适的阈值，对小波系数进行阈值处理，从而达到去除噪声而保留有用信号的目的.

2.3 算法描述

一方面，图像经过小波变换后，边缘和噪声是高频分量，图像的主要信息都集中在低频部分.由于各项异性扩散对噪声较敏感，而低频分量中的噪声较少，因此在低频部分进行各项异性扩散可以减小噪声对其的影响.另一方面，由参考文献［17］可知，小波收缩方法收敛速度快，而各项异性扩散需要较多次迭代，计算量较大.由于经过多尺度小波分解后，低频部分的大小远小于原图像.因此，本文在降噪部分采用混合去噪算法，即在高频分量部分采用小波收缩软阈值去噪，在低频采用基于差分曲率的各项异性扩散.算法的具体步骤如下：

1）低剂量投影数据是按式（1）的方法加在理想投影数据上，然后按式（2）进行MLEM算法重建.

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

3）对高频系数进行软阈值处理，去除噪声；对低频系数用基于差分曲率各项异性扩散算法进行图像降噪处理.

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

5）对各项异性扩散处理后的脉冲噪声进行中值滤波处理：由文献［18］可知，低剂量CT重建图像的噪声还表现为一些脉冲噪声，各项异性扩散降噪技术对重建的图像进行降噪后，仅可以平滑图像的小梯度区域，而相对于周围区域的大梯度区域则保持不变，这些大梯度可能是边缘，也可能是图像的峰值噪声.而中值滤波器只会对大噪声峰值产生的大梯度起作用，边缘产生的大梯度将不会受到影响.因此在低剂量CT重建时，低噪声可以由各项异性扩散平滑，而大噪声等脉冲噪声则由中值滤波器所消除.中值滤波公式为fn+1i，j=Median（fn+1i，j，w），其中w是中值算子的窗口.

6）重复步骤1）～5）直至得到最终的重建图像.

3 实验结果与讨论

为了验证算法的有效性，本文采用模拟人体骨骼组织和软组织的Hot-Cold模型作为实验模型进行低剂量CT图像重建的仿真，本文与BIMART（Block Iterative-Multiply algebraic reconstruction technique），BI-MLEM（Block Iterative-Maximum Likelihood Expe-ctation Maximization）［19］进行了比较，且为了进一步说明基于小波收缩和差分曲率的各项异性扩散算法在MLEM的低剂量CT重建算法的有效性，本文将传统P-M、方差应用到MLEM算法中，并和本文方法进行了比较.体模的大小为128 mm×128 mm，按照式（1）的方法向理想投影数据加入期望与方差的非线性关系的高斯噪声来仿真低剂量的投影数据，其中fi=200，T=1 200；实现本文算法的计算机操作系统为32位Microsoft Windows XP Professional 2002，处理器为英特尔奔腾双核E5300＠2.60 GHz，内存为2 G.计算结果如图1所示.

图1 各种算法重建图像的对比 Fig.1 Comparison of reconstructed image processed by various algorithms

由图1（b）可以明显看出，MLEM重建算法的重建结果图中含有较多噪声，即该算法不能有效地解决低剂量重建图像受噪声污染严重的问题.图1（c）和图1（d）与图1（b）相比可以看出：算法BI-MART和BI-EMML比MLEM算法重建后的图像包含更少的噪声，处理效果较好.将图1（g）的结果与图1（b）、（c）、（d）、（e）和（f）的处理结果相比可以看出：本文算法包含的噪声更少，图像更清晰，重建结果明显优于其它4种算法.综上所述，本文算法可以有效地解决低剂量重建图像的噪声问题，且可以在光滑去噪的同时较好地保持图像的纹理和边缘信息.

图2 肩部仿真模型重建 Fig.2 Comparison diagram of algorithm results

图2为肩部仿真模型的重建结果.将图2（e）与图2（a）、（b）、（c）和（d）相比，可以清晰地看出：（e）中包含的噪声较少，图像较平滑、清晰，图像质量较高，即本文算法处理后的图像明显优于其它算法处理后的图像.因此，本文算法可以有效地解决低剂量重建图像的噪声问题，且可以在光滑去噪的同时较好地保持图像的纹理和边缘信息.

从上述分析可得，本文提出的算法在去噪能力与保持边缘和细节能力方面，都明显优于其他重建方法.为了更好地证明算法的有效性，本文除了采用主观观察方法之外，对重建的结果也进行了客观、定量的分析.采用归一化均方误差、均方绝对误差、归一化均方距离、信噪比等评价参数对其进行评价.其定义分别为

1）归一化均方误差（Normalized Mean Squared Error，NMSE）

2）均方绝对误差（Mean Absolute Error，MAE）

3）归一化均方距离（Normalized Mean Square Distance，NMSD）

4）信噪比（Signal to Noise Ratio，SNR）

在公式（11）～（13）中，Fi和qi分别表示重建图像与原始图像的第i个像素的灰度值；Mi和mi分别表示重建图像与原始图像的均值；M和N分别为图像的行数和列数.

本文提出的新算法与各种算法比较的客观评价结果如表1所示.

表1 各种算法的客观评价 Tab.1 The objective evaluation of various algorithms

由表1可以明显看出，本文算法的归一化均方误差、均方绝对误差、归一化均方距离比其他几种方法都要小，说明了重建图像与原始图像的偏差和误差都较小，相似度较高.信噪比比其他的几种算法都大，说明重建图像中含有的噪声最少.由以上分析可知，本文算法得到的重建图像更接近原始理想图像.因此结合图2和表1可知，无论在视觉方面还是在定量评价方面，本算法在CT重建中是可行且有效的.

为了更好地比较4种算法的降噪效果，本文给出了肩部仿真模型的以上重建图像与原始理想图像的侧面轮廓线的比较图，选取第65行的灰度值进行比较，如图3所示.

图3 各种算法重建图像的截面图对比 Fig.3 Comparison of profiles of reconstructed image processed by various algorithms

从图3中可以清晰地看出，本文算法重建出的图像与理想图像的吻合度最高，更接近理想图像，即与其他几种算法相比，该重建图的噪声波动最小，去除噪声以及保留图像的细节信息的效果最好.

4 结 论

本文提出了一种基于小波变换和差分曲率各项异性扩散的MLEM的低剂量CT重建算法.该算法是采用基本的MLEM算法对低剂量投影数据进行重建，然后对重建后的图像在小波域进行降噪处理.通过与其他重建方法以及不同降噪方法的仿真实验相比较，无论是从主观的视觉效果还是从客观的质量评价，均表明本文提出的算法能够有效地去除噪声并保持图像的细节信息.算法可以应用于医学CT领域中，能够在降低辐射剂量的同时，有效地重建出符合要求的图像.因此，算法在医学成像领域具有较好的应用前景.

［1］唐瑜，崔学英，张权，等.基于偏微分方程的低剂量CT投影降噪算法［J］.计算机应用与软件，2014，31（3）：179-183.Tang Yu，Cui Xueying，Zhang Quan，et al.Denoising algorithm for low-dose CT projection based on partial differential equation［J］.Computer Applications and Software，2014，31（3）：179-183.（in Chinese）

［2］高志凌，刘祎，桂志国.基于模糊数学的低剂量CT投影域降噪算法［J］.测试技术学报，2011，25（6）：477-482.Gao Zhiling，Liu Yi，Gui Zhiguo.Noise reduction based on fuzzy theory for low-does CT sinograms［J］.Journal of Test and Measurement Technology，2011，25（6）：477-482.（in Chinese）

［3］刘祎，张权，桂志国.基于模糊熵的低剂量CT投影降噪算法研究［J］.电子与信息学报，2013，35（6）：1421-1427.Liu Yi，Zhang Quan，Gui Zhiguo.Noise reduction for low-dose CT sinogram based on fuzzy entropy［J］.Journal of Electronics ffInformation Technology，2013，35（6）：1421-1427.（in Chinese）

［4］Wang J，Li T F，Lu H B，et al.Penalized weighted least-squares approach to sinogram noise reduction and image reconstruction for low-dose X-ray computed tomography［J］.IEEE Transactions on medical imaging，2006，25（10）：1272-1283.

［5］Wang Jing，Lu Hongbing.Multiscale penalized weighted least-squares sonogram restoration for low-dose Xray computed tomography［J］.IEEE Transactions on Biomedical Engineering，2008，55（3）：1022-1032.

［6］李凯旋，黄静，马建华，等.低剂量CT重建中的双边滤波权值优化新方法［J］.电路与系统学报，2012，17（2）：94-99.Li Kaixuan，Huang Jing，Ma Jianhua，et al.Bilateral filtering weights optimization in low-dose computed tomography image reconstruction［J］.Journal of Circuits and Systems，2012，17（2）：94-99.（in Chinese）

［7］Zhang Q，Gui Z G，Chen Y，et al.Bayesian sinogram smoothing with an anisotropic diffusion weighted prior for low-dose X-ray computed tomography［J］.Optik-International Journal for Light and Electron Optics，2013，124（17）：2811-2816.

［8］高志凌.基于投影域降噪的低剂量CT重建算法研究［D］.太原：中北大学，2012.

［9］Rust G F，Aurich V，Reiser M.Noise dose reduction and image improvements in screening virtual colonoscopy with tube currents of 20 m As with nonlinear Gaussian filter chains［C］.Medical Imaging 2002 Conference.New York：IEEE，2002：186-197.

［10］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.

［11］Li T，Li X，Wang J，et al.Nonlinear sinogram smoothing for low-dose X-ray CT［J］.IEEE Transactions on Nuclear Science，2004，51（5）：2505-2513.

［12］张权，付学敬，桂志国，等.基于小波和各向异性扩散的PET图像MLEM重建算法［J］.计算机应用与软件，2013，30（11）：50-53.Zhang Quan，Fu Xuejing，Gui Zhiguo，et al.MLEMreconstrution algorithm for PET image based on wavelet and anisotropic diffusion［J］.Computer Applications and Software，2013，30（11）：50-53.（in Chinese）

［13］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.

［14］Chen Q，Montesinos P，Sun Q，et al.Adaptive total variation denoising based on difference curvature［J］.Image and Vision Computing，2010，28（3）：298-306.

［15］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.

［16］王知强.基于小波收缩与非线性扩散的去噪算法［J］.计算机工程，2011，37（7）：249-252.Wang Zhiqiang.Denoising algorithm based on wavelet shrinkage and nonlinear diffusion［J］.Computer Engineering，2011，37（7）：249-252.（in Chinese）

［17］Juan Liu，Moulin，P.Information-theoretic analysis of interscale and intrascale dependencies between image wavelet coefficients［J］.IEEE Transactions Image Processing，2001，10（11）：1647-1658.

［18］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.

［19］Byrne C L.Block-iterative methods for image reconstruction from projections［J］.IEEE Transactions on image processing，1996，5（5）：792-794.