This article presents the algorithms developed by the Core Imaging Library (CIL) developer team for the Helsinki Tomography Challenge 2022. The challenge focused on reconstructing 2D phantom shapes from limited-angle computed tomography (CT) data. The CIL team designed and implemented five reconstruction methods using CIL (https://ccpi.ac.uk/cil/), an open-source Python package for tomographic imaging. The CIL team adopted a model-based reconstruction strategy, unique to this challenge, with all other teams relying on deep-learning techniques. The iterative algorithms implemented with CIL showcased exceptional performance, with one algorithm securing third place in the competition. The best-performing algorithm employed careful CT data pre-processing and an optimization problem with single-sided directional total variation regularization combined with isotropic total variation and tailored lower and upper bounds. The reconstructions and segmentations achieved high quality for data with angular ranges down to 50 degrees, and in some cases acceptable performance even at 40 and 30 degrees. This study highlights the effectiveness of model-based approaches in limited-angle tomography and emphasizes the importance of proper algorithmic design leveraging on available prior knowledge to overcome data limitations. Finally, this study highlights the flexibility of CIL for prototyping and a comparison of different optimization methods.
Citation: |
Figure 4. Top Left: Absorption vs. path length for training dataset E. This shows the polynomial fit to pre-processed data and the estimated monochromatic correction. Top Right: Line profiles through the FDK reconstructions of dataset E with and without beam-hardening correction. Bottom left: FDK reconstruction of dataset A after beam-hardening correction has been applied. Bottom right: Line profiles through the FDK reconstruction of dataset A with and without beam-hardening correction
Figure 6. Reconstruction with FDK of the training sample A with 50$ ^\circ $ data. Left: reconstruction with linear gray color map in the range [-0.097, 0.104]. Center: segmentation. Right: Difference between segmentation and ground truth segmentation provided by the organizers. Correctly classified pixels are white, true background pixels misclassified as acrylic are red, and true acrylic pixels misclassified as background are blue
Figure 7. Reconstruction with isotropic TV of sample A with 50$ ^\circ $ data. Left: reconstruction displayed with linear gray color map in the range [-0.009, 0.064]. Center: segmentation. Right: Difference between segmentation and ground truth segmentation. Correctly classified pixels are white, true background pixels misclassified as acrylic are red, and true acrylic pixels misclasssified as background are blue
Figure 8. Fitting of the circle with the iterative method described in section 4.4, for the training sample B. In this example 4 iterations are required. In each plot we display the following: the fitted circumference, in red; the edges used by the algorithm in the current step to calculate the circle center and radius, in blue; the reconstruction of the complete data set; and the circle, whose radius is 4 pixels less than the fitted one, which determines which of the edges will be removed at the next iteration
Figure 9. Reconstruction with isotropic TV and disk constraint of sample A with 50$ ^\circ $ data. Left: reconstruction displayed with linear gray color map in the range [0.0, 0.0409]. Center: segmentation. Right: Difference between segmentation and ground truth segmentation. Correctly classified pixels are white, true background pixels misclassified as acrylic are red, and true acrylic pixels misclassified as background are blue
Figure 10. From left to right: isotropic TV reconstruction with a green arrow showing the direction of the center X-ray beam in the fan in the middle of the angular range for the limited angle dataset (i.e., the central direction), rotated isotropic TV reconstruction such that the y axis is parallel to the central direction (in green), isotropic plus anisotropic TV reconstruction as submitted to the challenge, isotropic plus anisotropic TV reconstruction at convergence (see section 4.8)
Figure 11. Reconstruction with isotropic and anisotropic TV with disk constraint of sample A with 50$ ^\circ $ data. The first row presents the results of the anisotropic TV algorithm as submitted to the challenge. The second row presents the results of the same algorithm at full convergence of the PDHG algorithm, resulting in an improvement of the score of about 4% with respect to the submitted algorithm (see section 4.8). The reconstructions are displayed with linear gray color map in the range [0.0, 0.0409]
Figure 13. Quantitative comparison of the five algorithms 24A–24E submitted by the CIL team (solid lines, with the best-performing algorithm 24B in black) and the best performing algorithm from each of the competing teams, 15A, 16B, 9, 8B, 17A, 14, 13 and 7 (dashed lines). Left: Average scores over the three phantoms at each of the seven levels from $ 90^\circ $ to $ 30^\circ $ data. Right: Scores transformed by $ -\log(1-\texttt{score}) $ to highlight differences between initial scores close to 1
Figure 14. Each column shows one of the evaluation samples, with ground truth on the top and output of algorithm 24B on the bottom. Superimposed as a green arrow is the central direction of the X-ray cone. Left: level 4 ($60^\circ)$, sample C. Center: level 4 ($60^\circ)$, sample B. Right: level 6 ($40^\circ$), sample C. If most edges are roughly not perpendicular to the central direction of the X-ray fan, these are captured by data, and score is high (Left: 0.96374). If most edges are roughly perpendicular, they are harder to reconstruct, and score is lower (center: 0.92751, right: 0.74646). The right column shows that only the structures parallel to the central direction of the X-ray can be recovered
[1] | J. Adler and O. Öktem, Learned primal-dual reconstruction, IEEE Transactions on Medical Imaging, 37 (2018), 1322-1332. doi: 10.1109/TMI.2018.2799231. |
[2] | A. H. Andersen and A. C. Kak, Simultaneous algebraic reconstruction technique (SART): A superior implementation of the ART algorithm, Ultrasonic Imaging, 6 (1984), 81-94. doi: 10.1016/0161-7346(84)90008-7. |
[3] | K. J. Batenburg and J. Sijbers, DART: A practical reconstruction algorithm for discrete tomography, IEEE Trans. Image Process., 20 (2011), 2542-2553. doi: 10.1109/TIP.2011.2131661. |
[4] | A. Beck and M. Teboulle, Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems, IEEE Trans. Image Process., 18 (2009), 2419-2434. doi: 10.1109/TIP.2009.2028250. |
[5] | M. Beister, D. Kolditz and W. A. Kalender, Iterative reconstruction methods in X-ray CT, Phys. Med., 28 (2012), 94-108, https://www.sciencedirect.com/science/article/pii/S112017971200004X. doi: 10.1016/j.ejmp.2012.01.003. |
[6] | J. Bian, J. H. Siewerdsen, X. Han, E. Y. Sidky, J. L. Prince, et al., Evaluation of sparseview reconstruction from flat-panel-detector cone-beam CT, Phys. Med. Biol., 55 (2010), 6575-6599. doi: 10.1088/0031-9155/55/22/001. |
[7] | A. Biguri, M. Dosanjh, S. Hancock and M. Soleimani, TIGRE: A MATLAB-GPU toolbox for CBCT image reconstruction, Biomed. Phys. Eng. Express., 2 (2016), 055010. doi: 10.1088/2057-1976/2/5/055010. |
[8] | C. Bouman and K. Sauer, A generalized Gaussian image model for edge-preserving MAP estimation, IEEE Trans. Image Process., 2 (1993), 296-310. doi: 10.1109/83.236536. |
[9] | R. A. Brooks and G. D. Chiro, Beam hardening in X-ray reconstructive tomography, Phys. Med. Biol., 21 (1976), 390. doi: 10.1088/0031-9155/21/3/004. |
[10] | R. Brown, C. Kolbitsch, C. Delplancke, E. Papoutsellis, J. Mayer, et al., Motion estimation and correction for simultaneous PET/MR using SIRF and CIL, Philos. Trans. A Math. Phys. Eng. Sci., 379 (2021), 20200208. doi: 10.1098/rsta.2020.0208. |
[11] | A. Chambolle and T. Pock, A first-order primal-dual algorithm for convex problems with applications to imaging, J. Math. Imaging Vis., 40 (2011), 120-145. doi: 10.1007/s10851-010-0251-1. |
[12] | T. F. Chan and S. Esedoglu, Aspects of total variation regularized l1 function approximation, SIAM Journal on Applied Mathematics, 65 (2005), 1817-1837. doi: 10.1137/040604297. |
[13] | G.-H. Chen, J. Tang and S. Leng, Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets, Med. Phys., 35 (2008), 660-663. doi: 10.1118/1.2836423. |
[14] | Z. Chen, X. Jin, L. Li and G. Wang, A limited-angle CT reconstruction method based on anisotropic TV minimization, Phys. Med. Biol., 58 (2013), 2119-2141. doi: 10.1088/0031-9155/58/7/2119. |
[15] | I. Coope, Circle fitting by linear and nonlinear least squares, J. Optim. Theory Appl., 76 (1993), 381-388. doi: 10.1007/BF00939613. |
[16] | M. J. Ehrhardt, P. Markiewicz, M. Liljeroth, A. Barnes, V. Kolehmainen, et al., PET reconstruction with an anatomical MRI prior using parallel level sets, IEEE Trans. Med. Imaging, 35 (2016), 2189-2199. doi: 10.1109/TMI.2016.2549601. |
[17] | S. Esedoḡlu and S. J. Osher, Decomposition of images by the anisotropic Rudin-Osher-Fatemi model, Commun. Pure Appl. Math., 57 (2004), 1609-1626. doi: 10.1002/cpa.20045. |
[18] | L. A. Feldkamp, L. C. Davis and J. W. Kress, Practical cone-beam algorithm, J. Opt. Soc. Am. A, 1 (1984), 612-619. doi: 10.1364/JOSAA.1.000612. |
[19] | J. Frikel and E. T. Quinto, Characterization and reduction of artifacts in limited angle tomography, Inverse Probl., 29 (2013), 125007, 21 pp. doi: 10.1088/0266-5611/29/12/125007. |
[20] | P. C. Hansen, J. Sauer Jørgensen and W. R. B. Lionheart (eds.), Computed Tomography: Algorithms, Insight, and Just Enough Theory, SIAM, Philadelphia, 2021. |
[21] | G. T. Herman, Correction for beam hardening in computed tomography, Phys. Med. Biol., 24 (1979), 81. doi: 10.1088/0031-9155/24/1/008. |
[22] | G. T. Herman and S. S. Trivedi, A comparative study of two postreconstruction beam hardening correction methods, IEEE Trans. Med. Imaging, 2 (1983), 128-135. doi: 10.1109/TMI.1983.4307626. |
[23] | S. Holman and P. Richardson, SPECT with a multi-bang assumption on attenuation, Inverse Probl., 36 (2020), 125005, 31 pp. doi: 10.1088/1361-6420/abab59. |
[24] | J. Hubbell and S. Seltzer, Tables of x-ray mass attenuation coefficients and mass energy-absorption coefficients, 2004. doi: 10.18434/T4D01F. |
[25] | J. S. Jørgensen, E. Ametova, G. Burca, G. Fardell, E. Papoutsellis et al., Core Imaging Library - part Ⅰ: A versatile Python framework for tomographic imaging, Philos Trans. A Math. Phys. Eng. Sci., 379 (2021), 20200192. doi: 10.1098/rsta.2020.0192. |
[26] | J. S. Jørgensen, E. Papoutsellis, L. Murgatroyd, G. Fardell and E. Pasca, CIL-HTC2022-Algo1, 2023. doi: 10.5281/zenodo.10212407. |
[27] | J. S. Jørgensen, E. Papoutsellis, L. Murgatroyd, G. Fardell and E. Pasca, CIL-HTC2022-Algo2, 2023. doi: 10.5281/zenodo.10212401. |
[28] | J. S. Jørgensen, E. Papoutsellis, L. Murgatroyd, G. Fardell and E. Pasca, CIL-HTC2022-Algo3, 2023. doi: 10.5281/zenodo.10212414. |
[29] | J. S. Jørgensen, E. Papoutsellis, L. Murgatroyd, G. Fardell and E. Pasca, CIL-HTC2022-Algo4, 2023. doi: 10.5281/zenodo.10212391. |
[30] | J. S. Jørgensen, E. Papoutsellis, L. Murgatroyd, G. Fardell and E. Pasca, CIL-HTC2022-Algo5, 2023. doi: 10.5281/zenodo.10212416. |
[31] | P. M. Joseph and R. D. Spital, A method for correcting bone induced artifacts in computed tomography scanners, J. Comput. Assist. Tomogr., 2 (1978), 100-108. doi: 10.1097/00004728-197801000-00017. |
[32] | 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, Phys. Med. Biol., 57 (2012), 7923-7956. doi: 10.1088/0031-9155/57/23/7923. |
[33] | Y. Lou, T. Zeng, S. Osher and J. Xin, A weighted difference of anisotropic and isotropic total variation model for image processing, SIAM J. Imaging Sci., 8 (2015), 1798-1823. doi: 10.1137/14098435X. |
[34] | A. Meaney, F. Silva de Moura, M. Juvonen and S. Siltanen, Helsinki Tomography Challenge 2022 (HTC2022) open tomographic dataset, 2023. doi: 10.5281/zenodo.8041800. |
[35] | S. Niu, Y. Gao, Z. Bian, J. Huang, W. Chen, G. Yu, Z. Liang and J. Ma, Sparse-view X-ray CT reconstruction via total generalized variation regularization, Phys. Med. Biol., 59 (2014), 2997-3017. doi: 10.1088/0031-9155/59/12/2997. |
[36] | N. Otsu, A threshold selection method from gray-level histograms, IEEE Trans. Syst. Man. Cybern., 9 (1979), 62-66. doi: 10.1109/TSMC.1979.4310076. |
[37] | E. Papoutsellis, E. Ametova, C. Delplancke, G. Fardell, J. S. Jørgensen, et al., Core Imaging Library - part Ⅱ: Multichannel reconstruction for dynamic and spectral tomography, Philos. Trans. A Math. Phys. Eng. Sci., 379 (2021), 20200193. doi: 10.1098/rsta.2020.0193. |
[38] | E. Pasca, J. S. Jørgensen, E. Papoutsellis, E. Ametova, G. Fardell, K. Thielemans, L. Murgatroyd, M. Duff and H. Robarts, Core Imaging Library (CIL), 2017-2023. doi: 10.5281/zenodo.4746198. |
[39] | E. Resmerita and R. S. Anderssen, Joint additive Kullback–Leibler residual minimization and regularization for linear inverse problems, Math. Methods. Appl. Sci., 30 (2007), 1527-1544. doi: 10.1002/mma.855. |
[40] | E. Y. Sidky, J. S. Jørgensen and X. Pan, Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm, Phys. Med. Biol., 57 (2012), 3065-3091. doi: 10.1088/0031-9155/57/10/3065. |
[41] | E. Y. Sidky, C.-M. Kao and X. Pan, Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT, J. Xray Sci. Technol., 119-139. |
[42] | E. Y. Sidky and X. Pan, Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization, Phys. Med. Biol., 53 (2008), 4777-4807. doi: 10.1088/0031-9155/53/17/021. |
[43] | E. Strekalovskiy and D. Cremers, Real-time minimization of the piecewise smooth Mumford-Shah functional, Lecture Notes in Computer Science, 8690 (2014), 127-141. doi: 10.1007/978-3-319-10605-2_9. |
[44] | Z. Tian, X. Jia, K. Yuan, T. Pan and S. B. Jiang, Low-dose CT reconstruction via edge-preserving total variation regularization, Phys. Med. Biol., 56 (2011), 5949-5967. doi: 10.1088/0031-9155/56/18/011. |
[45] | W. van Aarle, W. J. Palenstijn, J. Cant, E. Janssens, F. Bleichrodt, et al., Fast and flexible X-ray tomography using the ASTRA toolbox, Opt. Express., 24 (2016), 25129-25147. doi: 10.1364/oe.24.025129. |
[46] | S. Van der Walt, J. L. Schönberger, J. Nunez-Iglesias, F. Boulogne, J. D. Warner, et al., Scikit-image: Image processing in Python, Peer. J., 2 (2014), e453. |
[47] | C. Wang, M. Tao, J. G. Nagy and Y. Lou, Limited-angle CT reconstruction via the L1/L2 minimization, SIAM J. Imaging Sci., 14 (2021), 749-777. doi: 10.1137/20M1341490. |
[48] | T. Wang, K. Nakamoto, H. Zhang and H. Liu, Reweighted anisotropic total variation minimization for limited-angle CT reconstruction, IEEE Trans. Nucl. Sci., 64 (2017), 2742-2760. doi: 10.1109/TNS.2017.2750199. |
[49] | Z. Zhang, B. Chen, D. Xia, E. Y. Sidky and X. Pan, Directional-TV algorithm for image reconstruction from limited-angular-range data, Med. Image Anal., 70 (2021), 102030. doi: 10.1016/j.media.2021.102030. |
The 5 training datasets provided. The top row shows the full 360° sinograms A-E. The bottom row shows the ground truth reconstructions A-E
Left: Histogram of the pixel transmission values before and after re-normalization. Right: Zoom around 1.0
The sinogram, image data of the FDK reconstruction and the system geometry for the original (top) and zero-padded (bottom) data
Top Left: Absorption vs. path length for training dataset E. This shows the polynomial fit to pre-processed data and the estimated monochromatic correction. Top Right: Line profiles through the FDK reconstructions of dataset E with and without beam-hardening correction. Bottom left: FDK reconstruction of dataset A after beam-hardening correction has been applied. Bottom right: Line profiles through the FDK reconstruction of dataset A with and without beam-hardening correction
Performance comparison of Otsu and multi-Otsu segmentation on
Reconstruction with FDK of the training sample A with 50
Reconstruction with isotropic TV of sample A with 50
Fitting of the circle with the iterative method described in section 4.4, for the training sample B. In this example 4 iterations are required. In each plot we display the following: the fitted circumference, in red; the edges used by the algorithm in the current step to calculate the circle center and radius, in blue; the reconstruction of the complete data set; and the circle, whose radius is 4 pixels less than the fitted one, which determines which of the edges will be removed at the next iteration
Reconstruction with isotropic TV and disk constraint of sample A with 50
From left to right: isotropic TV reconstruction with a green arrow showing the direction of the center X-ray beam in the fan in the middle of the angular range for the limited angle dataset (i.e., the central direction), rotated isotropic TV reconstruction such that the y axis is parallel to the central direction (in green), isotropic plus anisotropic TV reconstruction as submitted to the challenge, isotropic plus anisotropic TV reconstruction at convergence (see section 4.8)
Reconstruction with isotropic and anisotropic TV with disk constraint of sample A with 50
Results of Algorithm 24B, compared to ground truth, for evaluation datasets
Quantitative comparison of the five algorithms 24A–24E submitted by the CIL team (solid lines, with the best-performing algorithm 24B in black) and the best performing algorithm from each of the competing teams, 15A, 16B, 9, 8B, 17A, 14, 13 and 7 (dashed lines). Left: Average scores over the three phantoms at each of the seven levels from
Each column shows one of the evaluation samples, with ground truth on the top and output of algorithm 24B on the bottom. Superimposed as a green arrow is the central direction of the X-ray cone. Left: level 4 (