2018, Vol. 39
计算机断层成像(computer tomography, CT)技术被公认为是20世纪下半叶最伟大的发明之一, 它给医学领域带来了革命性的变化, 还成功地应用于许多其他工程领域.对于CT成像系统而言，CT重建算法至关重要, CT图像重建算法的优劣直接影响着图像的质量与重建速度, 迄今为止已经出现了多种重建算法, 例如:FDK重建算法、BPF重建算法、ART算法等.这些算法各有优点，适用于不同条件下的成像, 但也都存在一些不足, 比如:FDK算法虽然重建速度快、占用内存少, 但锥角增大时伪影严重且反投影时间长；BPF算法虽然重建图像精确且可降低剂量, 但实现过程运算量大[1]；ART算法虽可重建角度缺失的投影数据，但重建时间较长[2].本文中引入了一种基于增广拉格朗日方法(augmented Lagrangian method)求解全变分正则化(total variation regularization)算法(ALMTVR), 此方法与其他一些经典重建算法相比有着重建速度快、重建图像质量高的优点.
1 增广拉格朗日方法的全变分正则化算法介绍一个用全变分正则化解决压缩感知问题的方法[3]为
|
(1) |
式中:u∈Rn或写成u∈Rs×t其中s×t=n, Diu∈R2表示图像u中第i个像素的离散梯度; A∈Rm×n(m < n)表示系统矩阵; b∈Rm是图像u经由系统矩阵A得到的测量值; ‖·‖可以代表 1范数(对应各项异性全变分), 也可以代表2范数(对应各项同性全变分).这种算法可以重建不同边界条件的图像u.
1.1 全变分最小化的增广拉格朗日算法在该ALMTVR算法中采用式(2)来等价替代式(1)中的全变分模型[4]:
|
(2) |
其相对应的增广拉格朗日函数为
|
(3) |
替换式(2)中的全变分模型后仍能保证全局收敛性定理收敛, 所以仍可以应用增广拉格朗日算法[5].根据式(3)中v′i与λ′的值须随着迭代次数不断变化, 将式(3)的真值记为u*和wi*, 则这两个因子的更新公式为
|
(4) |
|
(5) |
如果将增广拉格朗日算法直接应用到式(1), 则可以得到相应的增广拉格朗日函数为
|
(6) |
如果在每个像素中引入一个松弛因子wi∈R2, 则可以将Diu从不可微分的项‖·‖中分离出来, 然后去掉两者之间的差, 可以将式(6)中第一项求和中的每一项分解成
|
(7) |
把式(7)代回到式(6)中可以得到与式(3)中相同的目标函数.
根据增广拉格朗日算法的框架需要在每次迭代中有效地解出LA(wi, u)以求解式(1)中的问题.
1.2 整体算法流程通过将增广拉格朗日算法与交替最小化算法结合可以有效优化式(1)中的全变分模型, ALMTVR具体算法如下:
步骤1 初始化所有i值下的vi0, βi0, λ0, μ0和起始点wi0, u0;
步骤2 当不满足外部停止条件时:
1) 令
2) 求解从
3) 应用式(4)与式(5)中的因子更新公式获得vik+1, λk+1；
4) 选择新的惩罚因子βik+1≥βik和μk+1≥μk，直到满足外部停止条件.
外部停止条件可由以下两种方法确定:式(2)中的最优性条件已经近似达到; 相对变化
这种算法框架灵活, 可应用于诸多领域中.
1.3 参数选择由于应用仿真数据重建已知原图像, 所以可以用均方误差(mean squared error, MSE)来度量参数选择的好坏, 但对于实际的CT扫描数据来说, 由于原图像未知则不能采用均方误差来度量参数选择的好坏, 所以这里引入一个无需已知原图像的度量参数，称为度量标准Q(Metric Q).度量标准Q是一种无参考型度量参数[6], 它通过沿着图像中边缘等几何特征计算图像的锐利度与对比度来度量图像的质量.
为了方便起见，在这里简单介绍一下度量标准Q的计算方法.对于一幅图像I:[1, N1]×[1, N2], 首先将图像I分成M个互斥的片段{Pi}i∈{0, …, M}, 则每个片段的梯度矩阵Gi的计算公式为
|
(8) |
之后用Gi的奇异值si, 1>si, 2来计算一致性：
|
(9) |
则单片段方向的度量标准Q为
|
(10) |
对于最终的度量标准Q, 只有一致性的值超过特定阈值的片段才会对其产生影响.令S⊂[1, M]代表通过阈值片段的下标的集合, 则图像I的度量标准Q为
|
(11) |
为了验证本文提出的ALMTVR算法在CT图像重建领域中的应用价值, 实验采用5组CT影像数据来比较ALMTVR算法与传统的ART算法在CT重建中的重建速度与精度, 其中3组数据为仿真CT影像数据, 2组为核桃和兔子的实际扫描数据.本次实验的数据处理均在Matlab中进行.
2.1 实验 2.1.1 仿真实验模拟数据使用的是FORBILD的头部、胸部模体[7]和Shepp-Logan模型.原图像分辨率为256×256, 投影数据是采用512通道同时经360°每隔1°进行一次数据采集得到的, 投影数据加入2%的白噪声进行图像重建.
2.1.2 实际实验核桃的实际扫描数据是用口腔CT原型系统进行采集, X射线的管电压与管电流分别设置为70 kV与2 mA.采用512通道同时经一周360°进行368次均匀采集, 一周数据采集共用时18 s, 重建得到的图像分辨率为256×256.
兔子的实际数据是采用锥束CT实验系统进行扫描, X射线管电压和管电流分别为85 kV, 2.7 mA.X射线探测器为1 024通道, 一周360°采集300个投影数据，重建图像的分辨率为256×256.
2.2 实验结果及分析 2.2.1 仿真实验结果及分析仿真实验结果如图 1所示, 各图的对比度均设置在1到2之间.对比图 1c与图 1d, 图 1g与图 1h, 图 1k与图 1l可以看到, 相比ART算法得到的图像, ALMTVR算法得到的图像的噪声点更少, 通过观察得出ALMTVR算法重建出的图像抗噪声干扰的能力更强.
|
图 1 仿真数据重建实验结果 Fig.1 Experiment results of simulation data (a)—Shepp-Logan原图像；(b)—Shepp-Logan原图像上加入2%白噪声后的图像；(c)—Shepp-Logan ALMTVR算法重建图像；(d)—Shepp-Logan ART算法重建图像；(e)—FORBILD头部的原图像；(f)—在FORBILD头部原图像上加入2%白噪声后的图像；(g)—FORBILD头部ALMTVR算法重建图像；(h)—FORBILD头部ART算法重建图像；(i)—FORBILD胸部的原图像；(j)—在FORBILD胸部原图像上加入2%白噪声后的图像；(k)—FORBILD胸部ALMTVR算法重建图像；(l)—FORBILD胸部ART算法重建图像. |
为了更加直观地对比两种重建算法的重建效果, 将图 1a~图 1d中线上像素的像素值在同一张图中绘成曲线, 如图 2a所示, 将图 1e, 图 1f, 图 1g, 图 1h中线上像素的像素值在同一张图中绘成曲线, 如图 2c所示, 将图 1i, 图 1j, 图 1k, 图 1l中线上像素的像素值在同一张图中绘成曲线, 如图 2e所示.观察上述各图可以看到, ALMTVR算法重建出的图像的像素值更稳定，与原图像的像素值更贴近, 而ART算法重建出图像的像素值波动更大，受噪声影响更大.
|
图 2 像素值曲线和MSE、Metric Q对比图 Fig.2 Pixel value curves and MSE, Metric Q contrast diagram (a)—Shepp-Logan的像素值曲线；(b)—Shepp-Logan的MSE，Metric Q对比图；(c)—FORBILD头部像素值曲线；(d)—FORBILD头部MSE，Metric Q对比图；(e)—FORBILD胸部像素值曲线；(f)—FORBILD胸部MSE，Metric Q对比图. |
ART算法与ALMTVR算法中都存在可调节参数来调节算法的重建效果[8], 本次实验中ART算法使用的主要调节参数为迭代参数k, ALMTVR算法的主要调节参数为惩罚因子μ.实验中采用MSE与Metric Q来度量图像的质量, Metric Q取得最大值时表明图像质量最佳, MSE取得最小值时表明图像质量为最佳[9].为确认Metric Q是否能用来度量CT重建得到图像的质量[10], 本次实验中将改变ART与ALMTVR两种算法的主要参数进行多次重建, 并计算得到多组MSE与Metric Q值, 绘成曲线图显示在同一张图片中, 各组仿真实验结果如图 2b, 图 2d, 图 2f所示.
观察上述各图可知, 对于ART算法, MSE取得最小值与Metric Q取得最大值的迭代参数k相同均在k=1处.对于ALMTVR算法, 三组仿真数据重建图像MSE取得最小值的惩罚因子分别是6, 7, 6, Metric Q取得最大值时均为5.对于ALMTVR算法来说, 为了进一步确认采用Metric Q和MSE进行重建参数选择的差别, 分别采用Metric Q取得最大值时与MSE取得最小值时的μ值重建各组图像, 之后计算用不同μ值重建出各组图像的MSE值, 如表 1所示.从表 1可以看出, 对于同一组仿真数据采用两个不同参数μ进行重建的结果很接近.因此, 表明Metric Q可以用来度量CT重建图像的质量, 并用来作为重建关键参数选取的依据.
|
表 1 对于同一组仿真数据μ取值不同时MSE对比 Table 1 Comparison of MSE values for the same set of simulation data with different μ |
第一组实验为核桃数据实验, 用ART与ALMTVR两种算法进行重建, 并不断改变ART算法的关键参数k与ALMTVR算法的关键参数μ进行多次重建, 得到多个Metric Q值并绘制成曲线, 如图 3所示.Metric Q值取最大时代表图像的重建质量最高, 则由图 3可得，ART算法的最佳参数为k=1, ALMTVR算法的最佳参数为μ=5.
|
图 3 核桃数据Metric Q变化曲线 Fig.3 Metric Q curves of walnut data |
第二组实验为兔子标本数据实验, 与第一组实验过程相似, 用ART与ALMTVR两种算法进行重建, 并不断改变ART算法的关键参数k与ALMTVR算法的关键参数μ进行多次重建, 得到多个Metric Q值并绘制成曲线, 如图 4所示.由图 4可得，ART算法的最佳参数为k=1, ALMTVR算法的最佳参数为μ=7.
|
图 4 兔子标本数据Metric Q变化曲线 Fig.4 Metric Q curves of rabbit data |
用两种算法的最佳参数分别进行重建，得到核桃实验数据和兔子标本数据的实验结果如图 5, 图 6所示.
|
图 5 核桃数据实验结果 Fig.5 Results of the nut data (a)—ART算法重建结果；(b)—ALMTVR算法重建结果；(c)—图 5a与图 5b中线上像素点的像素值曲线. |
图 5a, 图 5b的对比度均设置在0到30之间, 对比图 5a与图 5b, 可以看到ALMTVR算法的重建效果比起传统的ART算法重建质量高, 图像中的噪声更少, 图 5c的像素值对比图中, ALMTVR算法的重建图像比ART算法得到的重建图像的像素值更稳定.
图 6a, 图 6b的对比度均设置在0到60之间, 对比图 6a与图 6b, 并观察图 6c的像素值对比图, 可以看到ALMTVR算法重建的图像质量要优于ART算法重建的图像质量.
重建时间如表 2所示.从表 2可以看出, 在重建相同投影数据时, ALMTVR算法重建时间要明显少于ART算法.
|
|
表 2 ALMTVR与ART算法的重建时间对比 Table 2 Comparison of reconstruction time for ALMTVR and ART |
综合分析，本文提出的ALMTVR算法与传统代数迭代ART算法重建结果相比较, 可以看出, 不论是在重建图像质量, 还是重建耗时上, ALMTVR算法均要优于ART算法.
3 结论本文提出的ALMTVR算法在CT图像重建领域的应用效果优秀.相对于传统的ART算法, ALMTVR算法重建出的图像噪声更少, 图像细节更加清晰, 并且显著提高了图像的重建速度, 有很高的实际应用价值.
| [1] |
伍绍佳, 陈皓, 廖丽, 等.
BPF重建算法的CUDA并行实现[J]. 集成技术, 2014(5): 61–68.
( Wu Shao-jia, Chen Hao, Liao Li, et al. CUDA parallel implementation of BPF reconstruction algorithm[J]. Integration Technology, 2014(5): 61–68. ) |
| [2] |
郑源彩, 潘晋孝, 孔慧华.
基于线性插值方法的ART重建算法研究[J]. 数学的实践与认识, 2013, 43(24): 80–84.
( Zheng Yuan-cai, Pan Jin-xiao, Kong Hui-hua. The ART reconstruction algorithm based on linear interpolation method research[J]. Mathematics Practice and Understanding, 2013, 43(24): 80–84. DOI:10.3969/j.issn.1000-0984.2013.24.012 ) |
| [3] |
Vandeghinste B, Goossens B, Holen R V, et al. Iterative CT reconstruction using shearlet-based regularization[C]// Society of Photo-Optical Instrumentation Engineers (SPIE). San Diego, 2012: 3305-3317.
|
| [4] |
Li C. An efficient algorithm for total variation regularization with applications to the single pixel camera and compressive sensing[D]. Houston: Rice University, 2010.
|
| [5] |
常小凯.
二次半定规划的增广拉格朗日算法[J]. 计算数学, 2014, 36(2): 133–142.
( Chang Xiao-kai. Augmented Lagrange algorithm for two semidefinite programming[J]. Computing Mathematics, 2014, 36(2): 133–142. ) |
| [6] |
Zhu X, Milanfar P.
Automatic parameter selection for denoising algorithms using a no-reference measure of image content[J]. IEEE Transactions on Image Processing, 2010, 19(12): 3116–3132.
DOI:10.1109/TIP.2010.2052820 |
| [7] |
Yu Z, Noo F, Dennerlein F, et al.
Simulation tools for two-dimensional experiments in X-ray computed tomography using the FORBILD head phantom[J]. Physics in Medicine & Biology, 2012, 57(13): 237–252.
|
| [8] |
Nilchian M, Vonesch C, Modregger P, et al.
Fast iterative reconstruction of differential phase contrast X-ray tomograms[J]. Optics Express, 2013, 21(5): 5511–5528.
DOI:10.1364/OE.21.005511 |
| [9] |
Wieczorek M, Frikel J, Vogel J, et al.
X-ray computed tomography using curvelet sparse regularization[J]. Medical Physics, 2015, 42(4): 1555–1565.
DOI:10.1118/1.4914368 |
| [10] |
Gao H, Qi X S, Gao Y, et al.
Megavoltage CT imaging quality improvement on TomoTherapy via tensor framelet[J]. Medical Physics, 2013, 40(8): 08191901–08191910.
|