Computed Tomography Reconstruction Algorithm Based on Relative Total Variation Minimization
-
摘要: 总变差(TV)最小算法是一种有效的CT图像重建算法,可以对稀疏或含噪投影数据进行高精度重建。然而,在某些情况下,TV算法会产生阶梯状伪影。在图像去噪领域,相对TV算法展现了优于TV算法的性能。鉴于此,将相对TV模型引入图像重建,提出相对TV最小优化模型,并在自适应梯度下降-投影到凸集(ASD-POCS)框架下设计对应的求解算法,以进一步提升重建精度。以Shepp-Logan、FORBILD及真实CT图像仿真模体进行重建实验,验证了该算法的正确性并评估了算法的稀疏重建能力和抗噪能力。实验结果表明,相对TV算法可以实现逆犯罪,可以对稀疏或含噪投影数据进行高精度重建,与TV算法相比,该算法可以取得更高的重建精度。
-
关键词:
- 计算机断层成像 /
- 相对总变差最小 /
- 稀疏重建 /
- ASD-POCS算法
Abstract: The total variation (TV) minimization algorithm is an effective CT image reconstruction algorithm that can reconstruct sparse or noisy projection data with high accuracy. However, in some cases, the TV algorithm produces stepped artifacts. The relative TV algorithm outperforms TV algorithm in the field of image denoising. In view of this, the relative TV model is introduced into image reconstruction, a relative TV minimum optimization model is proposed, and the corresponding solution algorithm is designed under the framework of adaptive gradient descent projection to the convex set (ASD-POCS) to further improve reconstruction accuracy. The reconstruction experiments were conducted with Shepp Logan, Forbild, and real CT image simulation models to verify the anti-crime ability of the algorithm and evaluate the sparse reconstruction and anti-noise abilities of the algorithm. The experimental findings reveal that the algorithm outperforms the TV method in terms of anti-crime capability and the ability to reconstruct sparse or noisy projection data with high precision. Compared with the TV algorithm, the algorithm can achieve higher reconstruction accuracy. -
图 11 不同模体在50个投影角度并于投影数据中加入方差为0.05的高斯白噪声条件下使用RTV和TV最小化重建算法重建结果的RMSE趋势曲线的比较
Figure 11. Comparison of RMSE trend curves of different phantom reconstruction results the using RTV and TV minimization reconstruction algorithms with 50 projection angles and adding Gaussian white noise with variance of 0.05 added to the projection data
1. Repeat main loop 13. $for{\text{ } }i = 1{\text{:} }{N_{{\rm{grad}}} }{\text{ do RTV - ASD loop} }$: 2. $ {{\boldsymbol{f}}}_{p}={{\boldsymbol{f}}}^{(n)}$ ${\boldsymbol{d} }{ {\boldsymbol{f} }^{(k)} } = {\nabla _f}\left\| { { {\boldsymbol{f} }^{(n)} } } \right\|_{{\rm{RTV}}}^{(k)}$ 3. ${\text{for } }i = 1{\text{ : } }{N_d}{\text{ do } }{\boldsymbol{f} } = {\boldsymbol{f} } + \beta { {\boldsymbol{A} }_{\boldsymbol{i} } }\displaystyle \frac{ { {g_i} - { {\boldsymbol{A} }_{\boldsymbol{i} } }{\boldsymbol{ \times f} } } }{ { { {\boldsymbol{A} }_{\boldsymbol{i} } }{\boldsymbol{ \times } }{ {\boldsymbol{A} }_{\boldsymbol{j} } } } }$ $ {\boldsymbol{d}}{{\boldsymbol{f}}^{(k)}} = {{{\boldsymbol{d}}{{\boldsymbol{f}}^{(k)}}} \mathord{\left/ {\vphantom {{{\boldsymbol{d}}{{\boldsymbol{f}}^{(k)}}} {{{\left\| {{\boldsymbol{d}}{{\boldsymbol{f}}^{(k)}}} \right\|}_2}}}} \right. } {{{\left\| {{\boldsymbol{d}}{{\boldsymbol{f}}^{(k)}}} \right\|}_2}}} $ 4. ${\text{for } }i = 1{\text{ : } }{N_i}{\text{ do if } }{{\boldsymbol{f}}_i} < 0{\text{ then } }{{\boldsymbol{f}}_i} = 0$ $ {{\boldsymbol{f}}^{(n)}} = {{\boldsymbol{f}}^{(n)}} - {d_\alpha }*{\boldsymbol{d}}{{\boldsymbol{f}}^{(k)}} $ 5. ${ {\boldsymbol{f} }_{{\rm{res}}} } = { {\boldsymbol{f} }^{(n)} }$ end for 6. ${{\boldsymbol{g}}^{(n)}} = {\boldsymbol{A}} \times {{\boldsymbol{f}}^{(n)}}$ 14. $ \nabla {{\boldsymbol{f}}_{{\text{RTV}}}} = {\left\| {{{\boldsymbol{f}}^{(n)}} - {{\boldsymbol{f}}_p}} \right\|_2} $ 7. ${d_p} = {\left\| {{{\boldsymbol{g}}^{(n)}}{\boldsymbol{ - }}{{\boldsymbol{g}}_{\boldsymbol{0}}}} \right\|_2}$ 15. ${\rm{if}}{\text{ } }\nabla { {\boldsymbol{f} }_{ {\text{RTV} } } } > {r_{\max} }*\nabla { {\boldsymbol{f} }_{ {\text{POCX} } } }{\text{ } }{\rm{and}}{\text{ } }{d_p} > \in {\text{ then} }$ 10. $\nabla {{\boldsymbol{f}}_{{\text{POCX}}}} = {\left\| {{{\boldsymbol{f}}^{(n)}} - {{\boldsymbol{f}}_p}} \right\|_2}$ ${d_\alpha } = {d_\alpha }*{\alpha _{{\rm{red}}} }$ 11. ${\text{if first iterarion , then }}{d_\alpha } = \alpha \cdot \nabla {{\boldsymbol{f}}_{{\text{POCX}}}}$ 16. $\beta = \beta \cdot {\beta _{{\rm{red}}} }$ 12. ${{\boldsymbol{f}}_p} = {{\boldsymbol{f}}^{(n)}}$ 17. ${\text{until}}\left\{ {{\text{stopping criteria}}} \right\}$ 18. ${\text{return } }\;{{\boldsymbol{f}}_{ {\rm{res} } } }$ 表 1 RTV算法和TV算法重建FORBILD模体的RMSE和SSIM比较
Table 1. RMSE and SSIM comparison of RTV and TV algorithms for FORBILD phantom reconstruction
项目 算法 投影个数 20 30 40 50 RMSE RTV 30.0×10-5 9.96×10-5 7.11×10-5 5.60×10-5 TV 5130×10-5 2530×10-5 620×10-5 590×10-5 SSIM RTV 0.9587 0.9792 0.9972 0.9992 TV 0.9293 0.9583 0.9943 0.9986 表 2 RTV算法和TV算法重建Shepp-Logan模体的RMSE和SSIM比较
Table 2. RMSE and SSIM comparison of the RTV and TV algorithms for Shepp-Logan phantom reconstruction
项目 算法 投影个数 20 30 40 50 RMSE RTV 8.010×10-5 4.791×10-5 3.349×10-5 2.586×10-5 TV 1610×10-5 530×10-5 160×10-5 80×10-5 SSIM RTV 0.9866 0.9948 0.9978 0.9995 TV 0.9812 0.9923 0.9956 0.9991 表 3 RTV算法和TV算法重建真实CT模体的RMSE和SSIM比较
Table 3. RMSE and SSIM comparison of the RTV and TV algorithms for CT phantom reconstruction
项目 算法 投影个数 20 30 40 50 RMSE RTV 0.0403 0.0294 0.0210 0.0120 TV 0.0467 0.0302 0.0218 0.0168 SSIM RTV 0.9213 0.9546 0.9743 0.9956 TV 0.8977 0.9387 0.9542 0.9842 表 5 RTV算法和TV算法在不同等级噪声下重建Shepp-Logan模体的RMSE和SSIM比较
Table 5. Comparison of RMSE and SSM of the RTV and TV algorithms for reconstructing Shepp-Logan phantom under different levels of noise
项目 算法 噪声等级 0.01 0.02 0.03 0.04 0.05 RMSE RTV 0.0015 0.0028 0.0035 0.0048 0.0055 TV 0.0065 0.0081 0.0102 0.0126 0.0148 SSIM RTV 0.9943 0.9823 0.9781 0.9653 0.9526 TV 0.9914 0.9814 0.9727 0.9594 0.9512 表 4 RTV算法和TV算法在不同等级噪声下重建FORBILD模体的RMSE和SSIM比较
Table 4. Comparison of RMSE and SSM of the RTV and TV algorithms for reconstructing FORBILD phantom under different levels of noise
项目 算法 噪声等级 0.01 0.02 0.03 0.04 0.05 RMSE RTV 0.0021 0.0029 0.0032 0.0043 0.0050 TV 0.0085 0.0107 0.0134 0.0162 0.0184 SSIM RTV 0.9923 0.9843 0.9812 0.9756 0.9712 TV 0.9836 0.9793 0.9754 0.9698 0.9652 表 6 RTV算法和TV算法在不同等级噪声下重建真实CT模体的RMSE和SSIM比较
Table 6. Comparison of RMSE and SSM of the RTV and TV algorithms for reconstructing CT phantom under different levels of noise
项目 算法 噪声等级 0.01 0.02 0.03 0.04 0.05 RMSE RTV 0.0188 0.0198 0.0205 0.0225 0.0233 TV 0.0194 0.0215 0.0233 0.0247 0.0286 SSIM RTV 0.9258 0.9042 0.8938 0.8814 0.8715 TV 0.9123 0.8926 0.8626 0.8612 0.8523 -
[1] PAN X, SIDKY E Y, VANNIER M. Why do commercial CT scanners still employ traditional, filtered back-projection for image reconstruction?[J]. Inverse Problems, 2008, 25(12): 1230009. [2] SIDKY E Y, PAN X. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[J]. Physics in Medicine & Biology, 2008, 53(17): 4777−4807. [3] YU L, YU Z, SIDKY E Y, et al. Region of interest reconstruction from truncated data in circular cone-beam CT[J]. IEEE Transactions on Medical Imaging, 2006, 25(7): 869−881. doi: 10.1109/TMI.2006.872329 [4] CANDES E, ROMBERG J, TAO T. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[J]. IEEE Transactions on Information Theory, 2006, 52(2): 489-509. [5] VOGEL C R. Iterative method for total variation denoising[J]. SIAM Journal on Scientific Computing, 1996, 17(1): 227−238. doi: 10.1137/0917016 [6] SIDKY E Y, EMIL Y, KAO C M, et al. Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT[J]. Journal of X-ray Science & Technology, 2006, 14: 119−139. [7] DONOHO D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52: 1289-1306. [8] LIU B, KATSEVICH A, YU H. Interior tomography with curvelet-based regularization[J]. Journal of X-ray Science and Technology, 2016, 25(1): 1−13. [9] SIDKY E Y, CHARTRAND R, BOONE J M, et al. Constrained TpV minimization for enhanced exploitation of gradient sparsity: Application to CT image reconstruction[J]. IEEE Journal of Translational Engineering in Health & Medicine, 2014, 2(6): 1−18. [10] CHAMBOLLE A, POCK T. An introduction to continuous optimization for imaging[J]. Acta Numerica, 2016, 25: 161-319. [11] MBOLLE A, POCK T. A first-order primal-dual algorithm for convex problems with applications to imaging[J]. Journal of Mathematical Imaging and Vision, 2011, 40(1): 120−145. doi: 10.1007/s10851-010-0251-1 [12] SIDKY E Y, JRGENSEN J H, PAN X. Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm[J]. Physics in Medicine & Biology, 2011, 57(10): 3065−3091. [13] XU L, YAN Q, XIA Y, et al. Structure extraction from texture via relative total variation[J]. ACM Transactions on Graphics, 2012, 31(6): 139. [14] RIGIE D S, RIVIèRE P L. Joint reconstruction of multi-channel, spectral CT data via constrained total nuclear variation minimization[J]. Physics in Medicine & Biology, 2015, 60(5): 1741. [15] QIAO Z, REDLER G, EPEL B, et al. A balanced total-variation-Chambolle-Pock algorithm for EPR imaging[J]. Journal of Magnetic Resonance, 2021, 328: 107009. doi: 10.1016/j.jmr.2021.107009 [16] 闫慧文, 乔志伟. 基于ASD-POCS框架的高阶TpV图像重建算法[J]. CT理论与应用研究, 2021,30(3): 279−289. DOI: 10.15953/j.1004-4140.2021.30.03.01.YAN H W, QIAO Z W. High order TPV image reconstruction algorithm based on ASD-POCS framework[J]. CT Theory and Applications, 2021, 30(3): 279−289. DOI: 10.15953/j.1004-4140.2021.30.03.01. (in Chinese). [17] XIN J, LIANG L, CHEN Z, et al. Anisotropic total variation minimization method for limited-angle CT reconstruction[C]//Proc. SPIE 8506, Developments in X-ray Tomography VIII, 85061C, 2012. https://doi.org/10.1117/12.930339. [18] YU H, WANG G. A soft-threshold filtering approach for reconstruction from a limited number of projections[J]. Physics in Medicine & Biology, 2010, 55(13): 3905−3916. [19] 乔志伟. 总变差约束的数据分离最小图像重建模型及其Chambolle-Pock求解算法[J]. 物理学报, 2018,67(19): 198701. DOI: 10.7498/aps.67.20180839.QIAO Z W. Total variation constrained data separation minimum image reconstruction model and its Chambolle-Pock algorithm[J]. Acta Physica Sinica, 2018, 67(19): 198701. DOI: 10.7498/aps.67.20180839. (in Chinese). [20] 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. [21] SIDDON R L. Fast calculation of the exact radiological path for a three-dimensional CT array[J]. Medical Physics, 1985, 12(2): 252−255. doi: 10.1118/1.595715 [22] GORDON R, BENDER R, HERMAN G T. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography[J]. Journal of Theoretical Biology, 1970, 29(3): 471−481. doi: 10.1016/0022-5193(70)90109-8 -