L1范数字典约束的感兴趣区域CT图像重建算法

L1范数字典约束的感兴趣区域CT图像重建算法
WU Junfeng;MOU Xuanqin
【摘要】针对现有全变分(TV)约束感兴趣区域(ROI)重建方法易产生块状伪影、细小结构丢失的问题,提出了一种L1范数字典稀疏约束的ROI低剂量CT医学图像重建算法.首先将ROI医学图像重建问题转化为最优化问题,以罚加权最小二乘函数为保真项,L1范数字典稀疏表示为约束项构建目标函数;然后将目标函数分解为图像更新和字典稀疏表示两个子优化问题,并交替求解上述两个子优化问题,实现ROI图像重建.胸腔模体仿真实验结果表明,在分别添加光子数为1×105、5×104和1×104泊松噪声投影情况下,与TV约束重建方法相比,图像结构相似度(SSIM)分别提高约0.103 5、0.113 1和0.125 8,峰值信噪比分别提高4.88、4.93和5.44 dB.山羊肺部实际CT扫描实验结果进一步证明,本文算法能够有效地去除块状伪影且较好的保留细小结构.
【期刊名称】《西安交通大学学报》
【年(卷),期】2019(053)002
【总页数】7页(P163-169)
【关键词】低剂量CT;感兴趣区域;字典学习;医学图像重建
【作者】WU Junfeng;MOU Xuanqin
【作者单位】
【正文语种】中文
【中图分类】TN911.73
X射线计算机断层成像(CT)已被广泛应用于医学临床诊断、工业无损检测、安全检查等诸多领域,是研究物体内部结构不可或缺的工具。

但是,在很多实际应用中,受辐射剂量、X射线束宽度和探测器尺寸的限制,投影数据往往会在探测器方向截断,该问题称之为感兴趣区域(ROI)重建[1-2]或内重建。

由于截断的投影数据不满足完备性条件,采用传统的滤波反投影算法进行重建,ROI区域内边缘会出现灰度值跌落和直流偏移现象,难以满足实际应用需求。

2008年,Courdurier等和Kudo等分别证明:若事先预知ROI内的一个子区域,则可利用凸集投影(POCS)算法对截断希尔伯特变换(THT)数据进行求逆,实现ROI精确重建[3-4]。

随后,Yu等和Jin等分别提出采用截断奇异值和连续奇异值分解算法对THT数据求逆[5-6]。

2017年,Liu等进一步提出一种推广的截断奇异值内重建方法[7]。

然而,上述ROI重建方法对其内部的子区域信息具有很强的依赖性,在许多重要的实际场合中,诸如稀有化石无损检测和灌注心脏CT成像,子区域信息很难获得。

为了解决前述问题,国内外研究人员提出了基于压缩感知的内重建方法。

Yu等首次证明若ROI为分段常值函数,则可通过全变分(Total Variation,TV)最小化模型实现其精确重建,并提出采用交替进行的有序子集联合代数迭代算法和梯度下降法分别最小化投影数据误差项和TV项,用以获得ROI重建图像[9],但该方法仅仅停留在理论研究和理想投影数据仿真实验阶段,而在实际的CT扫描过程中,投影数据通常含有噪声。

Xu等提出了一种TV约束的统计迭代ROI重建方法(SIRTV),并提出迭代软阈值算法极小化目标函数,初步解决了低剂量投影数据的ROI重建问题[10]。

然而,当噪声较大时,基于TV约束的重建结果常常存在图像边缘不清晰、表达细小图像结构能力差、块状伪影等问题。

字典是从大量的图像块样本中训练得到的一组过完备基,基于字典的稀疏表示比TV
稀疏变换能更好地适应不同图像的结构特征,从而在抑制噪声的同时能更多地保留图像细小结构11-13]。

为此,本文提出一种基于L1范数字典稀疏约束的统计内重建框架,并推导交替最小化算法求解模型,实现ROI图像重建。

最后,利用胸腔模体仿真低剂量投影数据和实际CT扫描山羊肺部投影数据,分别对本文所提重建算法的有效性进行验证和讨论。

1 基础知识
1.1 统计迭代重建
假设μ表示一幅二维待重建CT图像,其向量表示形式为μ=[μ1,…,μj,…μJ]T。

根据CT成像的物理过程,CT重建的数学模型可表示为
Aμ=P
(1)
式中:A=[ai,j]为I行、J列的系统矩阵,其中元素ai,j为加权因子,表示第j个像素对第i条射线线积分投影的贡献,1≤i≤I;P=[p1,…,pi,…pI]T为线积分投影数据或简称为投影数据。

统计迭代重建是一种有效的低剂量CT重建方法,该方法将CT图像重建问题转化为以待重建图像μ为变量的极小化问题[14]
(2)
式中:wi表示统计加权项,且wi=e-[Aμ]i;∀i=1,2,…,I,β为正则化参数;R(μ)为正则化项,表示待重建图像μ的先验信息。

1.2 稀疏表示和字典学习
设x∈RN×1表示大小为像素的图像块,D=[d1,d2,…,dM]∈RN×M(N≪M)表示字典,则根据图像稀疏表示理论,图像块x可以用该字典线性表示为[15]
(3)
式中:α=[α1,α2,…,αM]为稀疏编码系数。

由于N≪M,式(3)是欠定方程,有无穷多解,因此可通过求解如下极小化问题来实现
(4)
式中:λ为正则化参数。

为了构造好的字典,通常会选择与待重建图像相同类型的高质量图像中来提取图像块,构成训练样本。

字典学习的过程就是寻找一个能够使给定的训练样本集合中每个图像块都能用该字典中的原子稀疏表示的字典,即为求解如下优化问题
(5)
式中:X=[x1,…,xS]∈RN×S为选定的图像块样本;α=[α1,…,αN]为其相应的稀疏表示向量。

2010年,Mairal等提出一种快速在线学习法求解优化问题式(5),该方法在基于当前字典对图像块进行稀疏表示和字典更新两个步骤之间交替进行,直到算法满足停止迭代准则[16]。

本文选择稀疏模型软件优化工具箱进行字典学习[17]。

2 L1范数字典约束的ROI重建算法
2.1 目标函数建立
基于字典学习的稀疏表示在保持细节和抑制噪声方面均有上佳的表现,因此基于L1范数字典稀疏约束的统计迭代重建框架可描述为
(6)
式中:Rs表示图像块提取算子;αs表示第s个图像块的稀疏编码系数,s=1,2,…,S表示图像块数;β>0、λ>0均为正则化参数。

2.2 极小化目标函数
由于直接求解式(6)比较困难,因此本文采用交替最小化方法分别优化变量α和μ,即先固定变量α,求解目标函数关于变量μ的极小化问题;然后固定变量μ,求解关于变量α的极小化问题,不断重复上述交替求解过程,直到满足某种迭代收敛准则。

(1)固定α,更新μ。

当固定α时,优化问题(6)变换为如下形式
(7)
优化问题式(7)采用分离抛物面替代方法[14]求解,计算公式如下
(8)
式中:k为迭代次数。

为了加快算法的收敛速度,本文采用有序子集策略[14,18]。

(2)固定μ,更新α。

当固定μ时,优化问题(6)转化为如下形式
(9)
优化问题(9)等价于如下约束优化问题
(10)
式中:ε为稀疏编码误差。

式(10)为最小绝对收缩与选择算子模型(LASSO),可采用坐标下降[19]和最小角回归方法(LARS)[15]求解。

2.3 算法流程
综上,本文算法的流程如下:
(1)选择高质量CT图像并提取大小为8×8像素的图像块构造训练样本;
(2)设置字典学习参数,采用在线学习[16]方法训练包含256个原子字典D;
(3)设置参数ε、β、K(最大迭代次数)和L(有序子集个数);
(4)初始化μ1=0,α1=0,k=1;
(5)当没有达到最大迭代次数时,根据式(8)更新μk,l至μk,l+1;
(6)根据式(10)及LARS-LASSO算法求解αk,l+1关于字典D的稀疏表示;
(7)μk+1,1=μk,L,αk+1,1=αk,L,输出重建结果μ。

3 实验结果与分析
3.1 人体胸腔仿真实验
仿真实验选择胸腔模体生成一组理想投影数据,体模大小为512×512像素,每个像素的实际大小是0.976 6 mm×0.976 6 mm。

等角扇束在360°圆形轨迹下均匀采集1 160个视角的投影数据,X射线源到旋转中心和探测器的距离分别是570 mm 和1 400 mm,探测器单元数目为200。

为了仿真低剂量投影数据,向上述理想投影数据分别添加光子数为1×105、5×104、和1×104的泊松噪声,获得3组低剂量投影数据。

胸腔模体和选定的ROI如图1所示,胸腔模体的显示窗是[-1
000,900]HU。

图1 胸腔模体
本文首先从胸腔模体中提取大小为8×8像素的图像块作为训练样本,然后剔除方差特别小的图像块并删除训练样本中每个图像块的直流,最后使用在线学习方法训练含有256个原子的字典,结果如图2所示。

图2 利用图1训练出来的字典
编码误差设定为ε=2×10-7,有序子集数L为40,在分别添加光子数为1×105、
5×104和1×104泊松噪声投影下,正则化参数β分别设定为2×10-3、4×10-3和
8×10-3。

为了本文验证算法的收敛性,图3给出了添加光子数为5×104泊松噪声
投影情况下目标函数随迭代次数的变化曲线。

由图3可以看出,本文算法经过43次迭代后,目标函数值基本保持不变,因此本文设定迭代次数K=50时算法停止迭代。

图3 添加光子数为5×104泊松噪声投影情况下目标函数随迭代次数的变化曲线
两种算法在不同噪声水平下的重建结果如图4所示,其显示窗为[-1 000,900]HU。

可以看出,当噪声逐渐增大时,SIRTV重建结果有明显的块状伪影,而本文算法能够消除块状伪影且较好地保留细小结构(详细位置见图中箭头所指)。

1×105个光子1×105个光子
5×104个光子5×104个光子
1×104个光子1×104个光子(a)SIRTV算法 (b)本文算法图4 不同噪声水平下SIRTV算法和本文算法的重建结果
为了定量评估所提算法的重建性能,采用两种经典的评价指标对重建结果进行评估。

两种指标分别为图像结构相似度指数[20](SSIM,用符号ISSIM表示)和峰值信噪比(PSNR,用符号RPSN表示)。

表1列出了不同算法重建结果的结构相似度指数和峰值信噪。

表1结果可以说明本文算法的重建质量高于SIRTV算法。

表1 两种算法重建图像的结构相似度指数和峰值信噪比光子数ISSIMSIRTV本文
算法RPSN/dBSIRTV本文算法1×1050.832 20.935 720.9825.865×1040.810 10.923 218.4623.391×1040.732 90.858 715.8721.31
3.2 山羊肺部实际CT数据
为了进一步验证本文算法在实际CT数据的重建性能,利用西门子Somatom Sensation 64排CT对一头麻醉山羊进行正常剂量(球管电压为100 kV,球管电流
为150 mA)和低剂量(球管电压为80 kV,球管电流为17 mA)扫描,X射线源扫描轨迹半径为57 cm,在360°的扫描范围内均匀采集1 160个视角投影数据,等角探测
器包含672个探测单元,视野半径为25.05 cm。

为了获得穿过ROI区域的实际投
影数据,首先利用文献[11]算法对2组山羊肺部全局数据进行重建,重建图像为768×768像素,实际大小为43.63 cm×43.63 cm,并将其重建图像作为原始模体;然后选择半径为61像素的区域作为ROI;最后只保留穿过ROI的投影数据,即可获得仅穿过ROI的截断投影数据。

重建结果和选定的ROI如图5所示,其显示窗为[-700,800]HU。

(a)正常剂量 (b)低剂量图5 实际数据的重建结果
为了获得高质量字典,从图4a的肺部区域提取8×8像素的图像块作为训练样本,并采用3.1节所述字典学习方法训练包含256个原子的字典。

编码误差设定为ε=6.25×10-5。

在1 160、580和290个视角情况下,有序子集数分别设定为40、20和10,迭代次数K分别设定为50、80和100,正常剂量实验正则化参数β分别设定为0.015、0.022和0.033,低剂量实验正则化参数β分别设定为:0.020、0.030和0.060。

图6示出正常剂量山羊肺部数据不同视角数下两种算法的重建结果,图7给出1 160个视角两种不同算法重建ROI的水平中心和垂直中心剖线。

图8给出低剂量山羊肺部数据重建结果,图9给出低剂量山羊肺部数据1 160个视角下两种不同算法重建ROI的水平中心和垂直中心剖线,可以看出本文算法的细节保持能力优于
1 160个视角 1 160个视角
580个视角 580个视角
290个视角 290个视角(a)SIRTV算法 (b)本文算法图6 正常剂量山羊数据不同视角数情况下SIRTV算法和本文算法的重建结果
(a)水平中心
(b)垂直中心图7 1 160个视角下正常剂量山羊数据的重建ROI剖线图
SIRTV算法。

4 结论
1 160个视角 1 160个视角
580个视角 580个视角
290个视角 290个视角(a)SIRTV算法 (b)本文算法图8 低剂量山羊数据不同视角数情况下SIRTV算法和本文算法的重建结果
(a)水平中心
(b)垂直中心图9 1 160个视角下低剂量山羊数据的重建ROI剖线图
本文提出了一种基于L1范数字典稀疏约束的低剂量ROI统计迭代重建算法,不同于传统的TV稀疏约束,所提重建算法利用了高质量CT图像构造字典训练样本,采用了字典学习算法训练得到能够充分反映图像结构的字典,可以有效地保持图像细节信息。

胸腔模体仿真实验和山羊肺部实际CT扫描实验表明,所提重建算法能够有效去除块状伪影,保留重建图像细小结构。

致谢感谢马萨诸塞大学卢维尔分校俞恒永博士为本文提供了山羊肺部实际CT扫描数据。

参考文献:
【相关文献】
[1] WARD J P, LEE M, YE J C, et al. Interior tomography using 1D generalized total variation: part I Mathematical foundation [J]. SIAM Journal of Imaging Science, 2015, 8(4): 226-247.
[2] WANG G, YU H Y. The meaning of interior tomography [J]. Physics in Medicine and Biology, 2013, 58(16): 161-186.
[3] COURDURIER M, NOO F, DEFRISE M, et al. Solving the interior problem of computed tomography using a prior knowledge [J]. Inverse Problems, 2008, 24(6): 065001.
[4] KUDO H, COURDURIER M, NOO F, et al. Tiny a priori knowledge solves the interior problem in computed tomography [J]. Physics in Medicine and Biology, 2008, 53(9): 2207-2225.
[5] YU H Y, YE Y B, WANG G. Interior reconstruction using the truncated hilbert transform via singular value decomposition [J]. Journal of X-Ray Science and Technology, 2008,
16(4): 243-251.
[6] JIN X, KATSEVICH A, YU H Y, et al. Interior tomography with continuous singular value decomposition [J]. IEEE Transactions on Medical Imaging, 2012, 31(11): 2108-2119. [7] LIU R, HE L, LUO Y, et al. Singular value decomposition-based 2D image reconstruction for computed tomography [J]. Journal of X-Ray Science and Technology, 2017, 25(1): 113-134.
[8] YU H Y, WANG G. Compressed sensing based interior tomography [J]. Physics in Medicine and Biology, 2009, 54(9): 2791-2805.
[9] YU H Y, YANG J S, JIANG M, et al. Supplemental analysis on compressed sensing based interior tomography [J]. Physics in Medicine and Biology, 2009, 54(18): 425-432.
[10] XU Q, MOU X Q, WANG G, et al. Statistical interior tomography [J]. IEEE Transactions on Medical Imaging, 2011, 30(5): 1116-1128.
[11] XU Q, YU H Y, MOU X Q, et al. Low-dose X-ray CT reconstruction via dictionary learning [J]. IEEE Transactions on Medical Imaging, 2012, 31(9): 1682-1697.
[12] 王爱齐, 徐坤, 宋爱民. 基于局部自相似的字典学习图像去噪方法 [J]. 大连交通大学学报, 2017, 38(4): 192-195.
WANG Aiqi, XU Kun, SONG Aimin. Image denoising by locally self-similarity dictionary learning [J]. Journal of Dalian Jiaotong University, 2017, 38(4): 192-195.
[13] 宋长明, 王赟. 融合低秩和稀疏表示的图像超分辨率重建算法 [J]. 西安交通大学学报, 2018,
52(7): 18-24.
SONG Changming, WANG Yun. Super resolution reconstruction algorithm combined low rank with sparse representation [J]. Journal of Xi’an Jiaotong University, 2018, 52(7): 18-24.
[14] 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. [15] ELAD M, AHARON M. Image denoising via sparse and redundant representations over learned dictionaries [J]. IEEE Transactions on Image Processing, 2006,
15(12): 3736-3745.
[16] MAIRAL J L, BACH F, PONCE J, et al. Online learning for matrix factorization and sparse coding [J]. Journal of Machine Learning Research, 2010, 11(1): 19-60.
[17] JULIEN M. Spams-matlab-v2.5 [EB/OL]. (2014-07-04)[2017-05-30].http:∥spams-devel.gforge.inria. fr/downloads.html.
[18] 毛宝林, 陈晓朝, 孝大宇, 等. 基于全变分最小化和快速一阶方法的低剂量CT有序子集图像重建[J]. 物理学报, 2014, 63(13): 138701.
MAO Baolin, CHEN Xiaozhao, XIAO Dayu, et al. Ordered subset image reconstruction
studied by means of total variation minimization and fast first-order method in low dose computed tomograhpy [J]. Acta Physica Sinica, 2014, 63(13): 138701.
[19] LANGE K. Coordinate descent algorithms for LASSO penalized regression [J]. The Annals Applied Statistics, 2008, 2(1): 224-244.
[20] WANG Z, BOVIK A C, SHEIKH H R, et al. Image quality assessment: from error visibility to structural similarity [J]. IEEE Transactions on Image Processing, 2004, 13(4): 600-612.。