Inner product regularized multi-energy X-ray tomography for material decomposition

Salla-Maaria Latva-Äijö; Filippo Zanetti; Ari-Pekka Honkanen; Simo Huotari; Jacek Gondzio; Matti Lassas; Samuli Siltanen

doi:10.3934/ammc.2024001

Article Contents

2024, Volume 2, Issue 1: 1-16. Doi: 10.3934/ammc.2024001

This issue Previous Article Next Article Modeling the dynamics of actin and myosin during bleb stabilization

Inner product regularized multi-energy X-ray tomography for material decomposition

1.
Department of Mathematics and Statistics, University of Helsinki, Finland
2.
School of Mathematics, University of Edinburgh, UK
3.
Department of Physics, University of Helsinki, Finland

^*Corresponding author: Salla-Maaria Latva-Äijö
^*Corresponding author: Salla-Maaria Latva-Äijö

Received: June 2023

Revised: February 2024

Early access: March 2024

Published: March 2024

The first author is supported by the Doctoral Programme in Mathematics and Statistics (DOMAST).

Abstract / Introduction Full Text(HTML) Figure(11) / Table(4) Related Papers Cited by

Abstract

Multi-energy X-ray tomography is studied for decomposing three materials using three X-ray energies and a classical energy-integrating detector. A regularization term is used, which includes the inner products between the material distribution functions, penalizing any overlap of different materials. An interior point method is used to solve the underlying quadratic optimization problem; a previously developed preconditioner is extended to the case with three materials, while its theoretical properties are analyzed for any number of materials. The strategy is tested on real data of a phantom embedded with Na$ _2 $SeO$ _3 $, Na$ _2 $SeO$ _4 $, and elemental selenium. These selenium-based materials exhibit K-edges suitable for investigating the proposed method. It is found that the two-dimensional distributions of selenium in different oxidation states can be mapped and distinguished from each other with the proposed algorithm. The results have applications in material science, chemistry, biology and medicine.

Keywords:

Mathematics Subject Classification: Primary: 58F15, 58F17; Secondary: 53C35.

Citation:

FullText(HTML)

Figure 1. Spectrum of eigenvalues of the normal equations of the IPM, with and without preconditioner, for each iteration in the case $ N = 32 $, $ k = 3 $, $ \alpha = 150 $, $ \beta = 120 $

Download: Full-size image PowerPoint slide

Figure 2. The figure shows a schematic drawing of the phantom. The cuboid PMMA (grey) was drilled with three holes, each of which was filled with a mixture of a selenium compound [elemental Se (red), Na$ _2 $SeO$ _3 $ (green) and Na$ _2 $SeO$ _4 $ (blue)] and starch. The holes were capped with tissue paper (very light grey)

Download: Full-size image PowerPoint slide

Figure 3. Here is a photograph of the tomography setup. The detector in the picture is a different TimePIX-based model than the one used in this work. There is no visible difference between the materials

Download: Full-size image PowerPoint slide

Figure 4. The schematic drawing of the XAS-CT setup. The polychromatic X-rays produced by the X-ray tube are monochromatized with the spherically bent crystal analyzer. The sample to be imaged is illuminated by the monochromatized beam by moving it away from the Rowland circle so that the defocused beam completely covers it. The beam transmitted through the sample is recorded with a position-sensitive detector.[10]

Download: Full-size image PowerPoint slide

Figure 5. Two random example projection images of the sample. These projection images were taken with energies 12.658 keV and 12.662 keV

Download: Full-size image PowerPoint slide

Figure 6. Example sinogram of the data. Row 80 (81 in Matlab syntax) was picked from all the projection images of the sinograms. Overall, there were three different sinograms corresponding to the energies 12.658 keV, 12.662 keV, and 12.685 keV. These values correspond the low, middle, and high energies in our measurement model. The sinogram depicted here is produced from 12.662 keV projection images

Download: Full-size image PowerPoint slide

Figure 7. The mass attenuation coefficient, $ \frac{\mu}{\rho} $, as a function of photon energy for elemental Selenium Se, Na$ _2 $O$ _3 $, and PMMA. In the plot, the K-edges of Se and Na$ _2 $O$ _3 $ are clearly visible, which makes the decomposing of the materials easier with the IP method. In the energy window that we use (12.54-12.80 keV), which is marked with a black box in the figure, the attenuation values of PMMA are smooth, which means that it appears similar in every energy image. So, we can subtract it out from all the different energy sinograms by subtracting low energy sinogram from all the others

Download: Full-size image PowerPoint slide

Figure 8. Beta demonstration. First row: $ \beta $ = 0, second row: $ \beta $ = 8000, third row: $ \beta $ = 16000, fourth row: $ \beta $ = 19500, last row: ground truth. Notice that when a larger value of $ \beta $ is used, the reconstructed distribution of each material is sharpened

Download: Full-size image PowerPoint slide

Figure 9. Beta demonstration with colors. This demonstration is a colored visualization of the effect of increasing the $ \beta $-term. Here each color corresponds to a material. The first column shows all materials in the same figure, and the three following columns show only one material in each figure. First row: $ \beta $ = 0, second row: $ \beta $ = 8000, third row: ground truth. Notice that when a larger value of $ \beta $ is used, the reconstructed distribution of each material is sharpened

Download: Full-size image PowerPoint slide

Figure 10. Results of simulated data with measured realistic coefficients. The first row is the simulated result for the different selenium samples. Here we have totally ignored the plastic container, which was present in the practical real-world measurements. The second row shows the perfectly separated materials, which we have used as ground truths in the error calculations

Download: Full-size image PowerPoint slide

Figure 11. Reconstruction results. Row (a) represents the FBP reconstructions, with the phantom matrix. Row (b) shows the result with the FBP method, when we have subtracted one measurement (taken below all the K-edges) from all sinograms before reconstruction. This subtraction removes the PMMA phantom matrix from the reconstructions. Row (c): IP reconstructions (also with subtracted phantom matrix)

Download: Full-size image PowerPoint slide

Table 1. Attenuation coefficients $ c_{ij} $ for energy $ i $ and material $ j $

	Material 1	Material 2	Material 3
Energy 1	22.73	8.56	3.51
Energy 2	5.95	12.32	10.88
Energy 3	7.81	3.51	27.77

| Show Table

DownLoad: CSV

Table 2. Error table for different $ \beta $ values. We can see that the $ L_2 $-errors (12) are getting smaller, structural similarity is increasing, and the HaarPSI index is also improving when we increase the $ \beta $ regularization. We cannot increase $ \beta $ beyond the value of $ \alpha $ because the problem may then become non-convex. Unfortunately, the $ \beta $ regularization does not reduce the random noise spots around the actual target. Indeed, the inner product term penalizes the overlapping materials in different images, but it is not constructed to reduce the noise in the resulting images

Mean errors:	L2	SSIM	haarPSI	Corr.pix %	$ \alpha $	$ \beta $
First row:	0.60	0.25	0.14	96	20000	0
Second row:	0.56	0.26	0.17	96.4	20000	8000
Third row:	0.55	0.39	0.19	96.3	20000	16000
Fourth row:	0.59	0.41	0.21	95.6	20000	19500
Ground truth:	0	1	1	100

| Show Table

DownLoad: CSV

Table 3. Error measures for the simulated data

Material	Method	L2	SSIM	HPSI
Se	IP	0.27	0.87	0.42
SeO$ _3 $	IP	0.25	0.81	0.47
SeO$ _4 $	IP	0.20	0.89	0.54

| Show Table

DownLoad: CSV

Table 4. The X-ray attenuation coefficients for the different materials at different energies. The coefficients were determined in the article of Honkanen et al. [10]

Energy [eV]	Se	SeO$ _3 $	SeO$ _4 $
12658	8.47	2.44	3.51
12662	7.99	11.86	10.89
12685	7.42	8.00	27.78

| Show Table

DownLoad: CSV

Related Papers

Cited by

References

[1]	K. Batenburg, S. Bals, J. Sijbers, C. Kübel, P. A. Midgley, J. C. Hernandez, U. Kaiser, E. R. Encina, E. A. Coronado and G. Van Tendeloo, 3D imaging of nanomaterials by discrete tomography, Ultramicroscopy, 109 (2009), 730-740. doi: 10.1016/j.ultramic.2009.01.009.
[2]	K. Batenburg and J. Sijbers, DART: A practical reconstruction algorithm for discrete tomography, IEEE Transactions on Image Processing, 20 (2011), 2542-2553. doi: 10.1109/TIP.2011.2131661.
[3]	J. Baumann, Z. Kiss, Sven Krimmel, A. Kuba, A. Nagy, L. Rodek, B. Schillinger and J. Stephan, Advances in Discrete Tomography and its Applications, Springer, Birkhäuser Boston, 2007.
[4]	M. Colombo and J. Gondzio, Further development of multiple centrality correctors for interior point methods, Computational Optimization and Applications, 41 (2008), 277-305. doi: 10.1007/s10589-007-9106-0.
[5]	J. Gondzio, Interior point methods 25 years later, European Journal of Operational Research, 218 (2012), 587-601. doi: 10.1016/j.ejor.2011.09.017.
[6]	J. Gondzio, M. Lassas, S. Latva-Äijö, S. Siltanen and F. Zanetti, Material-separating regularizer for multi-energy x-ray tomography, Inverse Problems, 38 (2022), Paper No. 025013, 26 pp. doi: 10.1088/1361-6420/ac4427.
[7]	P. V. Granton, S. I. Pollmann, N. L. Ford, M. Drangova and D. W. Holdsworth, Implementation of dual-and triple-energy cone-beam micro-CT for postreconstruction material decomposition, Medical Physics, Wiley Online Library, 35 (2008), 5030-5042. doi: 10.1118/1.2987668.
[8]	G. Herman and K. Attila, Advances in Discrete Tomography and its Applications, Springer Science & Business Media, Birkhäuser Boston, 2008.
[9]	G. Herman and K. Attila, Discrete Tomography: Foundations, Algorithms, and Applications, Springer Science & Business Media, Birkhäuser Boston, 2012.
[10]	A. Honkanen and S. Huotari, Monochromatic computed tomography using laboratory-scale setup, Scientific Reports, 13 (2023), 363. doi: 10.1038/s41598-023-27409-6.
[11]	Y. Long and J. A Fessler, Multi-material decomposition using statistical image reconstruction for spectral CT, IEEE Transactions on Medical Imaging, 33 (2014), 1614-1626. doi: 10.1109/TMI.2014.2320284.
[12]	MATLAB, The MathWorks Inc. MATLAB version: 9.13.0 (R2022b), Natick, Massachusetts, United States, 2022. Available from: https://www.mathworks.com.
[13]	R. Reisenhofer, S. Bosse, G. Kutyniok and T. Wiegand, A Haar wavelet-based perceptual similarity index for image quality assessment, Signal Processing: Image Communication, 61 (2018), 33-43. doi: 10.1016/j.image.2017.11.001.
[14]	E. B. Saloman, J. H. Hubbell and J. H. Scofield, X-ray attenuation cross sections for energies 100 eV to 100 keV and elements Z = 1 to Z = 92, Atomic Data and Nuclear Data Tables, Elsevier, 38 (1988), 1-196. doi: 10.1016/0304-4076(88)90023-1.
[15]	R. Solem, T. Dreier, I. Gonçalves and M. Bech, Material decomposition in low-energy micro-CT using a dual-threshold photon counting x-ray detector, Frontiers in Physics, Frontiers Media SA, 9 (2021), 673843. doi: 10.3389/fphy.2021.673843.
[16]	Z. Wang, A. Bovik, H. Sheikh and E. Simoncelli, Image quality assessment: From error visibility to structural similarity, IEEE Transactions on Image Processing, 13 (2004), 600-612. doi: 10.1109/TIP.2003.819861.
[17]	S. J. Wright, Primal-Dual Interior-Point Methods, SIAM, 1997. doi: 10.1137/1.9781611971453.
[18]	W. Wu, H. Yu, P. Chen, F. Luo, F. Liu, Q. Wang, Y. Zhu, Y. Zhang, J. Feng and H. Yu, Dictionary learning based image-domain material decomposition for spectral CT, Physics in Medicine & Biology, IOP Publishing, 65 (2020), 245006. doi: 10.1088/1361-6560/aba7ce.
[19]	F. Zanetti and J. Gondzio, A new stopping criterion for Krylov solvers applied in interior point methods, SIAM J. Sci. Comput., 45 (2023), A703-A728. doi: 10.1137/22M1490041.