基于小波和非局部的全变差中值先验重建算法

张 芳，张 权，崔学英，董婵婵，刘 祎，孙未雅，白云蛟，桂志国，2＋

（1.中北大学 电子测试技术国家重点实验室，山西 太原030051；2.中北大学 仪器科学与动态测试教育部重点实验室，山西 太原030051）

0 引 言

低剂量CT［1］投影数据噪声模型的数据特点已经被广泛研究，其中Wang等通过对多个体模进行反复实验和分析，最终得出低剂量CT 投影数据经过对数变换后的均值和方差之间呈非线性递增关系，近似服从非平稳高斯分布的结论。针对投影数据统计特性的分析，由于最大似然期望最大算法 （MLEM）在重建过程中考虑了观测数据的统计特性，且具有非负性、全局收敛性和计数保持的特点，重建图像的效果得到了较好的改善，故而被广泛的使用。然而在实际中，当迭代次数达到一定后，随着迭代次数的增多重建图像的质量会出现棋盘效应，从而造成图像失真。基于Bayesian理论的最大后验 （maximum a posterior，MAP）法可以有效解决上述问题［2］，该方法既考虑了低剂量CT投影数据的统计特性，又通过对先验分布加入先验信息，使进行多次迭代后仍可以很好地抑制噪声且克服MLEM 重建算法收敛慢的缺点。MAP 重建的思想是在传统的MLEM 算法的基础上加入先验约束，从而达到抑制噪声、平滑图像、增强图像边缘的目的［3］。然而，基于传统的贝叶斯方法提供的先验信息是有限的，往往会使低剂量重建图像出现阶梯状伪影和过平滑现象［4］。本文将从此问题出发对低剂量CT 重建进行进一步的研究。

由于利用图像稀疏的先验知识能够在低剂量的条件下很好的重建图像，因此在不完备CT 投影数据的重建过程中对得到的梯度图像求最小化的l1范数，即运用全变差可以改善低剂量CT 重建图像的质量，该方法虽然提高了重建图像的抗噪声性能，但重建图像依然会出现图像边缘不清晰且存在块状伪影的问题。小波收缩由于其很好的时频局部特性，可以在去除噪声的同时保持图像的边缘和细节信息［5］，非局部也有很好的保持边缘细节的能力，本文受其启发，提出一种基于小波收缩和非局部的TV 中值先验重建算法。该算法在每次迭代过程中，用基于TV 的MP重建算法对低剂量投影数据直接进行重建，然后再通过小波变换用小波收缩和非局部对重建后的结果图进行处理。实验结果表明，本文算法明显改善了图像的质量，同时图像的信噪比也得到了很好的提高。

1 中值先验重建算法

Alenius等在1997年提出MRP （median root prior）算法。该算法具有边缘保持效果，其基本思想为把中值滤波器应用于重建迭代过程中，从而使得图像的像素值逐渐接近其邻域像素值的中值［5，6］。但是由于该算法不是真正意义上的MAP算法，只是一种经验公式。2003年，Hsiao等在MRP的基础上，提出一种具有优异边缘保持能力的真正MAP算法。

Hsicao定义的中值先验的目标函数如下式

式中：Φ（y｜f）——对数型的似然函数，R（f，m）——一种新的先验分布的目标函数。y 和f 分别表示观测数据向量和图像向量，m 是与f 具有相同维数的辅助向量。

先验分布的目标函数如下

2 基于小波收缩和非局部的TV 中值先验重建算法

2.1 基于TV的MP重建算法

由于MAP方法引入了图像的先验分布信息，将重建结果约束在正则空间，故改善了重建图像的质量，提高了图像的信噪比，在一定程度上保持了图像的细节信息［7］。但是基于传统的贝叶斯方法只可以提供有限局部先验信息，因此重建图像会过平滑，而且会出现阶梯状边缘的伪影。

由于利用图像稀疏的先验知识能够在欠采样的条件下很好的重建图像。而将图像进行离散梯度变换可得到稀疏图像，因此在低剂量CT 投影数据的重建过程中对得到的梯度图像求最小化的l1范数可以改善低剂量CT 重建图像的质量。而梯度的l1范数即为所谓的全变差。受此启发，本文将全变差［8－10］引入到MAP的目标函数中，从而提高图像的抗噪声性能，基于TV 的MP重建算法的目标函数为

式中：全变差 （total variation，TV）的连续形式为

对于N×M 大小的图像，其离散形式为

对目标函数进行求解可得

则图像第j个元素fj的TV 偏导数为

2.2 基于小波收缩和非局部的TV中值先验重建算法

尽管基于TV 的MP 重建算法取得了很好的效果，但是应用TV 最小化这一先验依然存在一些问题，对于低剂量的CT 投影数据，重建图像会出现边缘不清晰现象且存在着块状伪影，为了解决此不足，本文在基于TV 的MP重建算法的基础上，提出了一种基于小波收缩和非局部的TV 中值先验重建算法，从而进一步提高了重建图像的质量。

小波阈值去噪又称小波收缩［11］，是在最小均方误差下达到近似最优时对小波系数进行统一处理。图像进行小波变换后，图像的大部分能量集中在低频部分，高频部分包含了图像的噪声和边缘信息，且只集中了少部分能量，在文献 ［12］中提出噪声对应的小波系数幅值在噪声处较小，在边缘处较大，因此采用小波收缩的方法处理图像可以得到令人满意的效果。和正交小波变换相比，平稳小波变换具有平移不变性从而更适合用于图像的处理［13，14］。因此本文选用平稳小波变换对图像进行处理。其中2j为相应的尺度。由此分析，本文对基于TV 的MP 重建算法的每次迭代中，再进行3级离散平稳小波分解，然后对集中较多噪声的高频部分进行小波收缩处理，对低频成分的处理选取可较好保持细节和边缘的方法，从而保证图像的细节信息。

Buades等在2005 年利用图像包含众多相似结构的特性，提出了一种非局部均值 （nonlocal means，NLM）的图像去噪算法。该算法的基本思想是在全局范围内搜索当前像素点所在图像块相应的相似块，然后进行加权平均达到去噪的目的。其具体的计算方法如下

上式是以像素k、像素j 为中心像素点的邻域内所有像素灰度值的加权Euclidean距离，它为像素i和像素j 的邻域的一种相似度测量，其中，r 表示离散化噪声图像，为某向量的模值；r（Vi）和r（Vj）为灰度值向量，即表示像素i和j 周围的局部子块的像素集合。两者之间的相似性决定像素i和j 之间的相似性。

受其启发，本文选取非局部作为小波变换后低频成分的处理方法，实验表明该算法可以在降噪的同时有效的保持图像的边缘和细节，改善图像的信噪比，提高图像的抗噪声性能。

通过上述的分析，本文具体重建算法如下：

（1）基于TV 的MP重建

（2）在基于TV 的MP 重建算法的每次迭代中进行小波收缩处理：对低频成分采用保持图像边缘和细节很好的非局部进行降噪

对高频成分采用软阈值收缩，本文选取的阈值为长度对数阈值。最后使用平稳小波逆变换进行重构，得到优化的图像。

（3）重复以上过程一定次数后，得到最终的重建图像。

3 实验结果与分析

3.1 重建图像比较

本文算法及所有比较算法的实验仿真环境为：Windows 7旗舰版32位SP1 （DirectX 11）操作系统，英特尔Celeron （赛扬）E3300＠2.50GHz双核处理器，2GB内存。编程工具使用MATLAB7.6.0 （R2008a）。本文选取大小为128mm×128mm 的Sheep－Logan头部剖面图模型作为实验对象，如图1 （a）所示。本文选取大小为 （128×128）×（128×128）的系统矩阵，本实验采用平行投影方式，且在180个角度中均匀选取128个投影方向，每个投影方向分布128对探测器。本文按照下式向理想投影数据加入期望和方差成非线性关系的高斯噪声来仿真低剂量CT 投影数据

其中，i＝1，2…N 表示探测器信道，N 表示信道总数，λi表示第i个探测器获得投影数据的均值，σi2表示第i个探测器获得投影数据的方差，ki表示第i 个探测器的参数，T 表示用于描述扫描系统校准过程的系统参数，对于给定的CT采集系统，ki与T 是给定的。本文参数选取：ki＝200，T ＝12000。为了验证本文算法的有效性，将本文算法与传统的MLEM，OS－PML－OSL （ordered subsets－penalized maximum likelihood－one step late），MRP （median root prior），和在每次重建迭代中使用基于方差的各项异性扩散降噪的算法进行了比较。

图1 各种算法的对比结果

其中图1 （a）为原始图像。图1 （b）采用传统的MLEM 重建算法。图1 （c）采用基于有序子集的OSPML－OSL重建算法。图1 （d）使用本文介绍的MRP 重建算法。为了把本文算法和降噪效果很好的各项异性扩散算法进行比较，图1 （e）在每次迭代中使用基于方差的各项异性扩散降噪算法对重建图像进行降噪。图1 （f）为本文提出的算法。各种算法中涉及到的各种参数以及迭代次数，均为反复实验后得到的最优值。从结果图中可得，传统的MLEM 重建图像的质量最差；OS－PML－OSL 得到的图像虽然可以达到降噪的效果，却也模糊了图像的边缘和细节；MRP算法利用中值先验对重建图像的噪声进行了一定程度的消除且可获得比较清晰的图像，但是图像中存在一些明显的块状阴影；基于方差的各项异性扩散降噪处理虽然可达到不错的效果，但是依然存在少量的块状阴影；本文算法不仅可以有效的解决低剂量CT 图像的噪声问题，而且较好的保持了图像的边缘和细节信息，从视觉上分析，重建效果优于其它几种算法，初步说明该算法的有效性。

3.2 重建精度比较

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

（1）归一化均方距离 （normalized mean square distance，NMSD）

（2）均方绝对误差 （mean absolute error，MAE）

（3）信噪比 （signal to noise ratio，SNR）

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

表1 各种算法的客观评价

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

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

4 结束语

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

本文提出了一种基于小波收缩和非局部的全变差中值先验重建算法。该算法先在中值先验MP 算法的基础上，引入降噪性能优异的TV 方法，对目标函数进行修订，形成基于TV 的MP 重建算法；该算法可以很好的改善图像的质量，但是依然会有一些块状伪影；小波收缩可以在去除噪声的同时较好的保持图像的细节信息，非局部也可以很好的保持图像的细节，结合这几种方法的优势，提出了本文的算法，即在基于TV 的MP重建算法的每次迭代中，在小波变换后的小波域再进行小波收缩和非局部降噪。实验结果表明，该算法改善了图像的质量，提高了图像的抗噪声性能，在主观效果和客观效果来看，均说明了该算法的可行性。

［1］Bian J G，Yang K，Boone J M，et al.Investigation of iterative image reconstruction in low－dose breast CT ［J］.Physics in Medicine and Biology，2014，59 （11）：2659－2685.

［2］Ma J H，Zhang H，Gao Y，et al.Iterative image reconstruction for cerebral perfusion CT using apre－contrast scan induced edge－preserving prior ［J］.Physics in Medicine and Biology，2012，57 （22）：7519－7542.

［3］HE Lingjun，PAN Jinxiao，KONG Huihua.Adaptive regularized MAP of CT image reconstruction method ［J］.Computer Engineering and Applications，2011，47 （28）：198－200 （in Chinese）. ［何玲君，潘晋孝，孔慧华.自适应正则MAP 的CT 图像重建方法研究 ［J］.计算机工程与应用，2011，47（28）：198－200.］

［4］LI Xiaohong，ZHANG Quan，LIU Yi，et al.High quality median prior image reconstruction algorithm based on wavelet shrinkage and forward－and－backward diffusion ［J］.Journal of Computer Applications，2012，32 （12）：3357－3360 （in Chinese）.［李晓红，张权，刘祎，等.基于小波收缩和正逆扩散结合的优质中值先验图像重建算法 ［J］.计算机应用，2012，32 （12）：3357－3360.］

［5］Hsiao I T，Huang H M，Lin K J.A fast OS－type Bayesian reconstruction with an edge－preserving median prior［J］.Journal of Instrumentation，2009，4 （7）：P07012.

［6］Karali E，Koutsouris D.Towards novel regularization approaches to PET image reconstruction ［J］.Journal of Biosciences and Medicines，2013，1 （2）：6－9.

［7］Tang J，Rahmim A.Bayesian PET image reconstruction incorporating anato－functional joint entropy ［J］.Physics in Medicine and Biology，2009，54 （23）：7063－7075.

［8］Liu X W，Huang L H.A new nonlocal total variation regularization algorithm for image denoising ［J］.Mathematics and Computers in Simulation，2014，97：224－233.

［9］Bhadauria H S，Dewal M L.Medical image denoising using adaptive fusion of curvelet transform and total variation ［J］.Computers ＆ Electrical Engineering，2013，39 （5）：1451－1460.

［10］Jovanovski O.Convergence bound in total variation for an image restoration model［J］.Statistics ＆ Probability Letters，2014，90：11－16.

［11］Lee K，Vidakovic B.Semi－supervised wavelet shrinkage［J］.Computational Statistics ＆ Data Analysis，2012，56 （6）：1681－1691.

［12］Fierro M，Ha H G，Ha Y H.Noise reduction based on partialreference，dual－tree complex wavelet transform shrinkage ［J］.IEEE Trans Image Process，2013，22 （5）：1859－1872.

［13］Deng J M，Yue H Z，Zhuo Z Z，et al.A stationary wavelet transform based approach to registration of planning CT and setup cone beam－CT images in radiotherapy ［J］.Journal of Medical Systems，2014，38 （5）：1－9.

［14］Abdullah A J.Denoising of an image using discrete stationary wavelet transform and various thresholding techniques ［J］.Journal of Signal and Information Processing，2013，4 （1）：33－41.