邓露珍1,冯鹏1,何鹏1,陈绵毅1,2,张伟1,魏彪1,米德伶1
(1.重庆大学光电技术及系统教育部重点实验室,重庆,400044,中国;
2.伦斯乐理工大学生物医学工程学院医学成像中心,特洛伊,12180,美国)
摘要:针对最小化图像全变差(TV)方法中方向性的不足,本文以X和Y方向梯度变换为基础,引入对角方向(p/4)的梯度变换,从而进一步提高图像表示的稀疏性,并运用最速下降法进行求解。实验结果表明,在投影数较少的条件下,本文算法重构结果在均方根误差(RMSE)与通用质量指数(UQI)方面均优于代数迭代重建(ART)与TV方法,说明本文算法更适合稀疏投影情况下的重建。
关键词:计算机断层成像技术(CT);少剂量;对角图像总变差(DTV);压缩感知
在计算机断层成像(Computed Tomography,CT)研究领域中,少投影或“稀疏”投影重建一直是业内关注的热点[1]。对于CT重建而言,“稀疏”投影重建本质上是一个通过数值方法求解欠定方程最优解的问题,受限于算法本身在收敛性,抗噪性以及鲁棒性等方面所存在的局限,卷积反投影(Filter Back Projection,FBP)[2]、代数重建算法(Algebraic Reconstruction Techniques,ART)[3]、联立代数重建算法(Simultaneous Algebraic Reconstruction Techniques,SART)[4]等经典CT重建算法所得到的图像往往无法令人满意。压缩感知理论(Compressive Sensing,CS)[5]则指出,若被重建图像能够进行稀疏表示,就能够在较大概率条件下利用较少的投影求解得出。而重建图像质量的好坏,与图像的稀疏表示方法密切相关。
基于CS的诸多CT图像重建算法中[6-12],Pan等于2007年所提出的图像全变差(Total Variation,TV)最小化方法因为简单易行,效果良好而为广大研究人员所青睐[6]。TV方法使用X和Y方向的梯度算子来稀疏表示图像,虽然实现简单,但梯度算子只有水平垂直两个方向,如此并未充分利用CT图像中丰富的边缘和细节等方向信息。换言之,X和Y方向的梯度算子并不是最优的图像表示方法,TV算法依然有提升的空间。
基于此,本文在TV图像重建迭代无明显变化时引入对角方向(p/4)的梯度变换,力图借助多方向信息使重建中的CT图像获得更为稀疏的表达。以此为基础,提出了一种基于对角TV(Diagonal TV,DTV )的CT图像重建算法。该算法首先运用TV重建CT图像,经过多次迭代无明显变化之后,对图像进行对角方向的梯度变换并运用最速下降法[13]进行求解,得到最终重构图像。
本文章节安排如下,第一节是引言,第二节简要介绍CT和压缩感知的基本原理、TV及DTV的基本原理,以及本文所提算法的基本步骤。第三节对本文算法进行重建仿真,并与ART、TV算法的重建结果进行对比分析和讨论。第四节对全文进行总结。
1 原理和方法
1.1 CT和压缩感知的基本原理
CT是通过对物体进行不同角度下的射线投影而获取物体横截面信息的成像技术,它能够显示出物体的内部特性[14]。CT的模型可表示为:
(1)
上式中A为测量矩阵或投影矩阵,f为目标图像,p为投影数据。传统的CT图像重建算法通常使用解析算法FBP或迭代算法ART及SART等来重建CT图像,但在“稀疏”投影条件下上述算法均难以重建高质量的CT图像。
2006年,Candes、Tao和Donoho等提出的压缩感知理论突破了传统的奈奎斯特采样定理的瓶颈,对于可压缩的(稀疏)的信号通过低于奈奎斯特频率进行数据采样后,仍能够精确恢复原始信号[5]。其基本原理如下式所示:
(2)
上式中,为稀疏变换,为其对应的反变换。p为f在A的线性测量结果,在CT中即为投
影数据。若目标图像f是稀疏的,则通过正交变换,可将f压缩为y,y较f维数更少,但能量更为集中。
1.2 TV及DTV的基本原理
根据压缩感知的基本原理,2007年,芝加哥大学Pan等人提出了使用梯度变换来稀疏表示图像,并利用图像的全变差最小化作为约束的TV重建算法,获得良好效果。其定义如下:
(3)
这里,定义为重建图像f的总变差, 为f的像素值。s和t分别代表行和列,如图1所示。
DTV定义如下:
(4)
1.3本文算法
为了利用图像多方向的信息,本文在使用TV方法迭代到一定程度图像无明显变化之后使用对角方向的梯度变换来对图像进行稀疏表示,同时使用最速下降法进行最优化求解来重建CT图像。该算法可表示为:
(5)
(6)
我们采用最速下降法来求解式(5)和式(6),有:
(7)
(8)
其中m和n分别为TV和DTV在最速下降法中的迭代索引,M为TV算法迭代的总次数,和分别为TV和DTV的梯度下降步长。ART算法重构后得到的图像可以作为TV迭代的初始图像f 0,TV经过M次迭代得到f M作为DTV迭代的初始图像。gTV 为TV的梯度下降方向,表示为:
其中是一个已知的正整数,在接下来的算数中我们取。
该算法通过内、外两个循环进行迭代计算,外层循环执行ART迭代重建算法,内层循环首先进行TV约束,迭代M次重建结果无明显变化之后再进行DTV约束。
2.实验结果及讨论
为便于比较,本文采用图2所示头部模型作为重建对象,并利用ART、TV和本文算法(简称为TV+DTV)进行对比实验。假设CT系统是一代扫描系统,样本大小为256×256,投影角度为30个,扫描范围为1°~360°,按照如下方式进行扫描,ART和TV的迭代次数为1000次,本文算法首先用TV方法迭代600次,再用DTV迭代400次,重建结果如图3所示。
由图3不难看出,30个投影角度下,ART重建结果最差,出现诸多伪影,边缘处也较为模糊。而TV与本文算法之间整体上没有太大差异,但边缘处本文算法较TV算法更为清晰,内部均匀性上TV算法则略有优势。为进一步比较算法在更少的投影数据下的重建情况,我们在20个投影角度下进行重建,迭代次数依然是ART和TV均1000次,本文算法首先用TV方法迭代800次,再用DTV迭代200次,重建结果如图4所示。显而易见,图4所示重建结果质量均较图3有所降低,这说明随着投影数据的减少,重建结果有所下降,但相同条件下,本文算法依然是最优的,不仅边缘细节保持良好,内部ROI区域灰度分布也很均匀。
同时,本文采用均方根误差(Root Mean Square Error,RMSE)及UQI(Universal Quality Index)[15]来对重构图像的质量进行评价。RMSE定义如下:
(11)
式中M、N为图像的长和宽,表示原始图像的像素值,表示要评价的图像的像素值。
Wang以及Bovic等提出了UQI模型,此模型被广泛采用,并用于衡量图像在相关度、亮度及对比度的失真等方面。UQI范围是[-1,1],当原始图像和目标图像完全一致的时候,UQI的值为1。计算公式如下:
(12)
其中是原始图像均值,为待评价图像均值,为原始图像方差,为待评价图像方差,另
(13)
表1给出了前述30个投影角度下(30 Projection)和20个投影角度下(20 Projection)两组实验重建结果的RMSE与UQI。从表1可以看出,在20个投影角度下所有方法重建图像RMSE都超过在30个投影角度下所重建图像,而UQI则均小于在30个投影角度下所重建图像,这与图3、图4所示相吻合,说明投影数据的减少对所有重建算法均有影响,重建质量都略有下降。无论投影数据如何,TV+DTV算法重构所得RMSE小于ART和TV算法,UQI则大于ART和TV,说明整体上本文算法重建效果最佳,重构后图像与原始图像更加吻合。需要指出的是,当迭代数较少时,对角TV相对于传统TV并无显著优势,且在抗噪性能上仍有待提高。
3.结论
图像重建算法的优劣决定CT重构图像质量的优劣。针对TV算法中方向性的不足,本文在传统的对X和Y方向进行梯度变换的TV图像重建迭代无明显变化时引入对角方向的梯度变换,并运用最速下降法进行求解,以得到最优的重建结果。实验结果表明,本文算法在整体视觉效果上与TV方法差异不甚明显,但重构图像的边缘细节均优于ART和TV;且在投影数据较少的情况下,本文算法在RMSE与UQI方面均显著优于ART与TV,说明本文算法更适合稀疏采样情况下的重建。
参考文献
[1] Wang G, Yu H, De Man B. An outlook on X-ray CT research and development[J]. Chinese Journal of Medical Instrumentation, 2008, 32(3): 157-169.
[2] Zhang Shun-li, Li Wei-bin, Tang Gao-feng. Study on Image Reconstruction Algorithm of Filtered Backprojection[J]. Journal of Xianyang Normal University, 2008, 23(4): 47-49.
[3] Gordon R, Bender R, Herman GT. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography[J]. Journal of theoretical Biology, 1970, 29(3): 471-481.
[4] Andersen H, Kak AC. Simultaneous algebraic reconstruction technique (SART): A superior implementation of the ART algorithm[J]. Ultrasonic imaging, 1984, 6(1): 81-94.
[5] Donoho DL. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306.
[6] Sidky EY, Kao C, Pan X. Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT[J]. Journal of X-Ray Science and Technology, 2006, 14: 119-139.
[7] E. Y. Sidky and X. C. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” Physics in medicine and biology, vol. 53, pp. 4777, 2008.
[8] Y. Liu, J. Ma, Y. Fan and Z. Liang, “Adaptive-weighted total variation minimization for sparse data toward low-dose x-ray computed tomography image reconstruction,” Physics in medicine and biology, vol. 57, pp. 7923, 2012.
[9] Y. Liu, Z. Liang, J. Ma, H. Lu, K. Wang, H. Zhang and W. Moore, “Total variation-stokes strategy for sparse-view x-ray CT image reconstruction,” IEEE Transactions on Medical Imaging, 2014
[10] M. Chen, D. Mi, P. He, L. Deng, and B. Wei, “A CT Reconstruction Algorithm Based on L 1/2 Regularization,” Computational and mathematical methods in medicine, 2014.
[11] M. Chang, L. Li, Z. Chen, Y. Xiao, L. Zhang, and G. Wang, “A few-view reweighted sparsity hunting (FRESH) method for CT image reconstruction”, Journal of X-ray science and technology, vol. 21, no. 2, pp. 161-176, 2013.
[12] L. Deng, P. Feng, M. Chen, P. He, Q. S. Vo and B. Wei, “A CT Reconstruction Algorithm Based on Non-Aliasing Contourlet Transform and Compressive Sensing,” Computational and Mathematical Methods in Medicine, 2014.
[13] 李化欣,潘晋孝. 最速下降法在图像重建中的应用[J]. 科技情报开发与经济. 2006, 16(3):155-156.
[14] 李新彩. 基于压缩感知的CT迭代图像重建技术应用研究[D]. 山东: 山东大学信号与信息处理, 2011: 1-60.
[15] Liu YR. Research on Objective Full-Reference Image Quality eva luation Method[J]. Nanjing: Computer Science & Technology College, 2010: 1-64.
评论信息