Open Access
Add to BibSonomyAbstract
A penalized maximum-likelihood estimation is proposed to perform hyperspectral (spatio-spectral) image reconstruction for X-ray fluorescence tomography. The approach minimizes a Poisson-based negative log-likelihood of the observed photon counts, and uses a penalty term that has the effect of encouraging local continuity of model parameter estimates in both spatial and spectral dimensions simultaneously. The performance of the reconstruction method is demonstrated with experimental data acquired from a seed of arabidopsis thaliana collected at the 13-ID-E microprobe beamline at the Advanced Photon Source. The resulting element distribution estimates with the proposed approach show significantly better reconstruction quality than the conventional analytical inversion approaches, and allows for a high data compression factor which can reduce data acquisition times remarkably. In particular, this technique provides the capability to tomographically reconstruct full energy dispersive spectra without compromising reconstruction artifacts that impact the interpretation of results.
© 2015 Optical Society of America
1. Introduction
Hard X-ray fluorescence (XRF) tomography has grown into a powerful non-destructive technique for probing trace metal distributions in samples without the need for physical sectioning that may disrupt trace element distributions and sample structure. The technique is becoming increasingly important in a growing number of research fields including the life sciences [1–3], the environmental and earth sciences [4–5], and the materials science [6]. A key advantage of XRF imaging is its high sensitivity and specificity for transition metals such as iron, copper, zinc, or other essential trace elements [7], its high resolution that allows imaging at scales from few microns [8] down to tens of nanometers [9–10], and for biological samples in particular, the capacity for imaging trace elements in whole, hydrated samples when cryogenic temperatures are used [11]. The results published to date have largely used conventional filtered back-projection reconstruction methods. While these have yielded very good results, the ultimate quality and precision of the inversion is dependent on both theoretical and experimental factors, such as dealing with finite sampling of data, the uncertainties in data, models, and perhaps more importantly, incorporation of prior knowledge about the nature of possible problem solutions [12].
In a typical XRF experiment, the specimen is illuminated by a focused monochromatic X-ray beam to stimulate emission of fluorescence X-rays from elements whose binding energies lie below the incident beam energy (see Fig. 1 for the schematic representation of XRF data collection). The fluorescent X-rays on the beam trajectory are considered to be emitted isotropically and the escaped photons are typically measured using solid state energy-dispersive X-ray detectors placed close to the specimen. Unlike full-field tomographic methods, the focused nature of the incident beam necessitates collecting image data through two-dimensional (2D) raster scanning of the sample. In tomographic mode, an energy-resolved dispersive spectra is acquired for each projection angle as a function of the horizontal, and in some cases also vertical, translation of the sample through the X-ray beam. An XRF tomographic experiment therefore consists of raster scans over multiple projection angles. Traditionally, the energy-dispersive data collected by each detector pixel is converted into separate elemental maps before performing a 3D tomographic reconstruction. This is usually performed by a curve fitting approach or a hyperspectral analysis [13–15]. These approaches, however, often fail for pixels with low photon counts. In this paper, we propose an alternative approach by reconstructing the spectrum at each spatial point from the complete spatio-spectral data, and to perform the 2D elemental mapping on the reconstructed hyperspectral slices. A detailed illustration of the proposed method is illustrated in Fig. 2. The advantage of this approach is that we can implement reconstruction algorithms imposing additional prior knowledge on the third, in this case energy, axis.
Fig. 1 A schematic representation of a typical XRF data collection system. The sample is illuminated by a focused monochromatic X-ray beam (for example using a Fresnel zone plate or elliptically figured X-ray reflective mirrors) to stimulate emission of fluorescence X-rays from elements whose binding energies lie below the incident beam energy. The fluorescent X-rays on beam trajectory are considered to be emitted isotropically and the escaped photons are typically measured using solid state energy-resolving X-ray detector placed close to the specimen. X-ray attenuation though the sample can be measured using a variety of different transmission detectors placed downstream of the sample (ion chambers, photodiodes, etc.)
Fig. 2 Schematic representations of the traditional and the proposed joint XFM elemental reconstruction approaches. The drawing on the left illustrates the traditional XFM data analysis, where each 1D energy-dispersive pixel data is converted into individual elemental maps (1) before performing a 3D tomographic reconstruction (2). The drawing on the right depicts a joint approach that uses the complete spatio-spectral data to perform reconstruction in a single step. In this way, the variations in both spatial (3) and energy (4) dimensions can be regularized in the reconstruction domain simultaneously, which in turn leads to a significantly reduced reconstruction uncertainty, and better interpretation of results.
To date, analytical inversion methods such as the filtered back projection (FBP) [16], or the regridding algorithm [17] have been most often used for the tomographic reconstruction of elemental distributions. A drawback of these methods is that the reconstructions often suffer from artifacts due to limited acquisition angles and poor signal-to-noise ratios, as is often the case of XRF [18]. Statistical approaches, on the other hand, not only allow for accurate modeling of the underlying physics of photon transport, and the corresponding photon detection statistics, but also allow regularization of the problem when the data is incomplete. Thus, they provide a versatile tool for dealing with low-photon counts and non-ideal data collection geometries [19–21]. In this paper, we present a general formalism of reconstructing hyperspectral images from the full spatiospectral XRM dataset. With this approach, rather than maximizing the data likelihood, Bayes’ theorem is used to maximize the posterior density which includes prior terms that enforce desired (e.g. smoothness, sharpness) and/or certain (e.g. non-negativity) properties of the solution [22]. This approach leads to objective functions that combine data misfit and prior terms, which can be optimized using various iterative [23] and global search methods [24]. In this paper, we focused on the feasibility of the hyperspectral inverse model and ignored self-absorption and scattering effects [25–26]. We present the performance of the reconstruction method with experimental data acquired from a seed of arabidopsis thaliana collected at the 13-ID-E microprobe beamline at the Advanced Photon Source.
2. Methods
2.1. Mathematical formulation of the problem
We cast the problem as a multidimensional inverse problem of reconstructing the full energy dispersive spectra for each voxel in the sample from measured data. We define xE ∈ ℝN as the virtual spectra per voxel, and yE ∈ ℝM as the spectrally resolved data. The corresponding problem is modeled as a discrete optimization problem in the form of:
where the objective function consists of the Poisson log-likelihood term ℒ(xE, y), and a penalty term ℛ(xE). The problem formulated in this way is commonly interpreted as a penalized maximum likelihood (PML) estimator [27], and is closely related with the maximum a posteriori estimator in Bayesian statistics [28]. This approach allows us to separately enforce regularizations of solutions in both the spatial and the spectral dimensions. A Poisson model is assumed to deal with low photon counts of XRF, and ℛ(xE) is used to regularize roughness in solution domain. The scalar β controls the amount of applied regularization, and is selected heuristically (typically from a range of 1 to 100) so it provides the desired image. The penalty term follows in the general form of a Gibbs prior [29] of, where V represents a set of voxels in a neighborhood of the ith voxel. For instance, in a two-dimensional (2D) space, it corresponds to the eight neighboring voxels in the spatial domain, and two neighboring voxels in the spectral domain (see Fig. 3 for details). The constants wij are nonnegative weights that satisfy wij=wji for all i and j, and ψ is a user-defined function to regulate the problem depending on the local statistics of model parameters. The weights are set inversely proportional to the distance between the neighboring voxels in a single dimension. For ψ, we used the function [30] of, which provides smooth solutions for large δ values, while favoring clustering or piecewise-constant solutions, for finite values of δ. This allows either a quadratic or a pseudo-linear model to describe parameter variations in their respective domains. A natural choice for the XRF inversion is to use larger δ values to impose smoothness in the spectral domain, and to use smaller δ for the spatial domain. A closed form expression for the solution of the optimization problem was obtained using the optimization transfer principle [31] (see, Appendix for details).Fig. 3 A three-dimensional ”data cube” (with 2D spatial and 1D spectral variations of model parameters), and the interactions among neighboring voxels in spatial and spectral dimensions are illustrated at the top. The regularization function for various δ values are plotted at the bottom.
2.2. Experimental data
The experimental dataset was collected at the GSECARS X-ray microprobe beamline, 13-ID-E, at the Advanced Photon Source, at Argonne National Laboratory. A single 350 μm diameter seed of arabidopsis thaliana was glued to a 100 μm diameter quartz fiber for analysis, mounted hanging down vertically from the rotational stage above the sample. An incident beam energy of 12 keV from a Si(111) double crystal monochromator was focused to about 2 × 2 μm with reflective rhodium-coated silicon mirrors in a Kirkpatrick-Baez geometry. A four-element silicon-drift diode detector array was used (Hitachi ME-4), coupled to a high-speed digital spectrometer system (XIA XMap), to measure the X-ray fluorescence signal. To collect the X-ray fluorescence data the seed was rotated in the focused beam through 363 degrees at an angular speed of 33.333 degrees per second, saving spectra accumulated through 0.5 degrees every 15 milliseconds. After each full rotation, the seed was stepped horizontally across the beam by 1 mm, and then rotated 363 in the opposite direction, generating spectral arrays consisting of 2048 energy bins, 727 angles, 501 horizontal steps, and 4 detector elements, which we will refer to here as our measured sinogram. The data were stored in HDF5 format [32].
2.3. Computational details
The proposed reconstruction algorithm is implemented on a high-performance data-intensive computing middleware [33] and performed image reconstruction at Argonne Leadership Computing Facility. Specifically, we ran our experiments on Mira, a 10-petaflops IBM Blue Gene/Q system. We used up to 896 nodes where each node consists of 16 physical cores (14336 cores in total) and 16 GiB memory. The running time of the total hyperspectral reconstruction takes about 6.77 seconds for 1 sPML iterations.
3. Results
3.1. Hyperspectral reconstruction results
The energy spectra were reconstructed using both the filtered-backprojection (FBP), and the spatio-spectral PML (sPML) method, and the corresponding results are presented in Fig. 4. The resolutions of the dataset were 1 μm and 15 eV, in spatial and energy dimensions, respectively. For the sPML reconstruction, δ was selected as 0.01 for the spatial dimensions. The penalty term with this set of parameters is expected to favor piecewise-constant solutions in the spatial coordinates, but still enforces some smoothness when the differences between neighboring voxels fall below a certain threshold, thus helping to further regulate the measurement noise and outliers. A larger δ of 0.1 was used to impose further smoothing on the reconstructed spectra. The locally obtained spectrum with sPML after 80 iterations is remarkably smooth compared to the one produced by FBP, which renders an accurate local elemental analysis. The sPML approach preserves many elemental fingerprints in the reconstructed spectrum otherwise are masked by the noise with FBP. The smaller values (< 10−1) are particularly hard to resolve at the pixel level due to the streaking artifacts that FBP produces. For example, the reconstructed values around 2.3 keV by FBP are all negative, and produce very inaccurate estimations compared to sPML. We used the main Kα 1 emission lines of iron (Fe), manganese (Mn) and zinc (Zn) elements to examine the reconstructed elemental maps. The energy windows for the Fe, Mn and Zn emissions were initially determined as 5.71 – 6.08 keV, 6.13 – 6.67 keV and 8.49 – 8.82 keV, respectively, and the corresponding overlaid images were obtained by integrating the highlighted regions of the spectrum. Streaking artifacts in the FBP reconstruction are visible where some localized elements are present at high levels of concentration. Such streaking artifacts propagate through the image reconstruction and degrade image quality, whereas the sPML is very robust to such streaking artifacts. A detailed comparison of individual elements is presented in Fig. 5. The corresponding cross-sections of the reconstructions are plotted.
Fig. 4 The absolute values of the reconstructed energy spectrum (in an arbitrary location inside the sample) using FBP (blue), absolute value of FBP (green) and sPML (red) are plotted on the top. The corresponding overlaid elemental maps of Mn, Fe, and Zn are shown on the bottom. The size of the seed is about 250 μm.
Fig. 5 Individual 2D elemental maps of Mn, Fe, and Zn and their cross-section plots for both FBP and sPML.
3.2. Effect of angular sampling
In order to reduce the data as well as experimental time, we have analyzed the effect of angular sampling of projections. The total number of projections for the complete dataset was 720, taken at 0.5 degree intervals from 0 to 360 degrees. We reduced the projections to simulate undersampled radial acquisitions with acceleration factors of 2, 4, an 8, that corresponds to 1, 2, and 4 degrees tilt intervals, respectively. The corresponding effect of undersampling on spatial resolution and stability of images for the cases of FBP and sPML is presented in Fig. 6. As an example, we presented a slice reconstruction at Fe-peak energy, which consists of many small details and Fe clusters. The tiny Fe clusters scattered around the large clusters (observable in the zoomed-in region) are distinctly resolvable with acceleration factors up-to 4 using sPML, whereas the aliasing artifacts which are typically enhanced with FBP due to undersampling suppress subtle details in reconstructed images, even with a small acceleration factor of 2. These artifacts gets worse with increasing acceleration factors for FBP, and makes the interpretation of results almost impossible. However, sPML is robust to noise and outliers even with high data compression ratios.
Fig. 6 The effect of angular sampling on reconstructed images at Fe-peak energy is presented for both reconstruction methods. The columns from right to left show the reconstructions obtained using a tilt step angle of 0.5, 1, 2, and 4 degrees, respectively. Bottom row shows the difference images between sPML and FBP. All the images are scaled to the same intensity range indicated by the color bar, and only the non-negative values are presented.
3.3. Relative error
Relative error (RE) was used to quantify the bias and errors induced by reconstruction methods. It was calculated as the following relative difference between the estimated and measured data,
where e represents the energy dependence. [xi] and [di] are, respectively the two-dimensional reconstructed image and the measured sinogram data with spatial indices represented by i. P is the number of projections. The relative error plot as a function of energy for different reconstruction methods (with 15 eV energy resolution) are presented in Fig. 7. The relative error of the proposed sPML algorithm is about 0.002 for the entire energy range, and is a significant indicator of reliable reconstructions. FBP, on the other hand, provides a relative error of almost two orders of magnitude higher than sPML. The main contributing factor for sPML having much better relative error is that, first, it doesn’t suffer from streaking artifacts as FBP does (mainly due to the non-negativity constraint); and second, the ability to suppress image noise by using appropriate regularization.4. Conclusions
In this paper, we introduced a general framework to perform hyperspectral image reconstruction for XRF tomography. This multidimensional inversion approach provides the capability to tomographically reconstruct the full energy dispersive spectra for each spatial voxel in the sample, and avoids conventional reconstruction artifacts that impact interpretation of results. It allows for a high data compression factor which can significantly reduce data acquisition times, and improve beamtime efficiency. We demonstrated the validity of the method with experimental data acquired from an arabidopsis seed sample. The estimated element distributions show good data fidelity and better image quality than the conventional analytical inversion approaches, and is robust to noise and outliers in measurements. The sPML algorithm will be made publicly available as part of the TomoPy software package [34]. We are currently extending this reconstruction approach to other multidimensional problems such as time-resolved tomography, or diffraction tomography.
Appendix
Optimization transfer refers to the methods which replace the original objective function, like given in expression 1, at each iteration with a surrogate function, which when maximized, is guaranteed to also increase the value of the original objective function. By choosing the surrogate functions appropriately, reductions in computation time and improvements in convergence can be realized. Here we only give the update equations of the algorithm, but for details one can refer to [31].
where γ(x) := ψ̇(x)/x (note that ψ(x) is given in equation 3). It is implicitly assumed that F(k) = 0 if β = 0. When β = 0, the iteration is equivalent to the ML expectation-maximization (ML-EM) solution.Acknowledgments
We thank Amanda Socha and Tracy Punshon (Dartmouth College) for sharing the arabidopsis thaliana data used in this paper. We also thank Stefan Vogt, Chris Jacobsen, Eugene Lavely, and Yi-San Lai for fruitful discussions and helpful comments during the course of the work. This research used resources of the U.S. Department of Energy (DOE) Office of Science User Facilities operated for the DOE Office of Science by Argonne National Laboratory under Contract No. DE-AC02-06CH11357. The GSECARS 13-ID-E beamline is supported by the National Science Foundation: Earth Sciences ( EAR-1128799), and Department of Energy: Geosciences ( DE-FG02-94ER14466).
References and links
1. T. Paunesku, S. Vogt, J. Maser, B. Lai, and G. Woloschak, “X-ray fluorescence microprobe imaging in biology and medicine,” J. Cell Biochem. 99, 1489–1502, (2006). [CrossRef] [PubMed]
2. C. J. Fahrni, “Biological applications of X-ray fluorescence microscopy: exploring the subcellular topography and speciation of transition metals,” Curr. Opin. Chem. Biol. 11, 121–127, (2007). [CrossRef] [PubMed]
3. S. Kim, T. Punshon, A. Lanzirotti, L. Li, J. Alonso, J. Ecker, J. Kaplan, and M. Guerinot, “Localization of iron in arabidopsis seed requires the vacuolar membrane transporter VIT1,” Science 314, 1295–1298, (2006). [CrossRef] [PubMed]
4. E. Lombi and J. Susini, “Synchrotron-based techniques for plant and soil science: opportunities, challenges and future perspectives,” Plant Soil 320, 1–35, (2009). [CrossRef]
5. E. Lombi, M. D. Jonge, E. Donner, C. G. Ryan, and D. Paterson, “Trends in hard X-ray fluorescence mapping: environmental applications in the age of fast detectors,” Anal. Bioanal. Chem. 400, 1637–1644, (2011). [CrossRef] [PubMed]
6. S. Bohic, A. Simionovici, X. Biquard, G. Martinez-Criado, and J. Susini, “Synchrotron X-ray microfluorescence and microspectroscopy: Application and perspectives in materials science,” Oil Gas Sci. Technol. 6, 979–993, (2005). [CrossRef]
7. M. D. Jonge and S. Vogt, “Hard X-ray fluorescence tomography – an emerging tool for structural visualization,” Curr. Opin. Struc. Biol. 20, 606–614, (2010). [CrossRef]
8. D. Bourassa, S. C. Gleber, S. Vogt, H. Yi, F. Will, H. Richter, C. H. Shin, and C. J. Fahrni, “3D imaging of transition metals in the zebrafish embryo by X-ray fluorescence microtomography,” Metallomics 9, 1648–1655, (2007).
9. H. Suhonen, F. Xu, L. Helfen, C. Ferrero, P. Vladimirov, and P. Cloetens, “X-ray phase contrast and fluorescence nanotomography for material studies,” Int. J. Mater. Res. 103, 179–183, (2012). [CrossRef]
10. H. M. Hertz, J. C. Larsson, U. Lundstrom, D. H. Larsson, and C. Vogt, “Laboratory x-ray fluorescence tomography for high-resolution nanoparticle bio-imaging,” Opt. Lett. 39, 2790–2793, (2014). [CrossRef] [PubMed]
11. S. Chen, J. Deng, Y. Yuan, C. Flachenecker, R. Mak, B. Hornberger, Q. Jin, D. Shu, B. Lai, J. Maser, C. Roehrig, T. Paunesku, S. C. Gleber, D. J. Vine, L. Finney, J. VonOsinski, M. Bolbat, I. Spink, Z. Chen, J. Steele, D. Trapp, J. Irwin, M. Feser, E. Snyder, K. Brister, C. Jacobsen, G. Woloschak, and S. Vogt, “The Bionanoprobe: Hard X-ray fluorescence nanoprobe with cryogenic capabilities,” J. Synchrotron Radiat. 21, 66–75, (2014). [CrossRef]
12. P. C. Hansen, “Discrete Inverse Problems: Insight and Algorithms,” : (SIAM, Philadelphia, PA2010).
13. A. J. Brown, “Spectral curve fitting for automatic hyperspectral data analysis,” IEEE Geosci Remote S 44, 1601 (2006).
14. A. J. Brown, B. Sutter, and S. Dunagan, “The MARTE VNIR imaging spectrometer experiment: design and analysis,” Astrobiology 8, 1001–1011, (2008). [CrossRef] [PubMed]
15. X. R. Wang, A. J. Brown, and B. Upcroft, “Applying incremental EM to Bayesian classifiers in the learning of hyperspectral remote sensing data,” Information Fusion Proc. 1, 606–613, (2005).
16. A. C. Kak and M. Slaney, “Principles of Computerized Tomographic Imaging,” : (SIAM, Philadelphia, PA2001).
17. B. A. Dowd, G. H. Campbell, R. B. Marr, V. Nagarkar, S. Tipnis, L. Axe, and D. P. Siddons, “Developments in synchrotron x-ray computed microtomography at the National Synchrotron Light Source,” Proc. SPIE 3772, 224–236, (1999). [CrossRef]
18. G-F. Rust and J. Weigelt, “X-Ray Fluorescent Computer Tomography with Synchrotron Radiation,” IEEE T. Nucl. Sci. 45, 75–88, (1998). [CrossRef]
19. J. Kaipo and E. Somersalo, “Statistical and Computational Inverse Problems,” : (Springer, New York, NY2004).
20. P. J. La Riviere and P. A. Vargas, “Monotonic penalized-likelihood image reconstruction for x-ray fluorescence computed tomography,” IEEE T. Med. Imaging 25, 1117–1129, (2006). [CrossRef]
21. P. J. La Riviere, P. A. Vargas, M. Newville, and S. Sutton, “Reduced-scan schemes for X-ray fluorescence computed tomography,” IEEE Trans. Nucl. Sci. 54, 1535–1542, (2007). [CrossRef]
22. J. Qi and R. M. Leahy, “Iterative reconstruction techniques in emission computed tomography,” Phys. Med. Biol. 51, R541 (2006). [CrossRef] [PubMed]
23. E. X. Miqueles and A. R. De Pierro, “Iterative reconstruction in X-ray fluorescence tomography based on Radon inversion,” IEEE T. Med. Imaging 30, 438–450, (2011). [CrossRef]
24. P. L. Riviere, P. Vargas, M. Rivers, and S. R. Sutton, “Penalized-likelihood image reconstruction for X-ray fluorescence computed tomography,” Opt. Eng. 45, 077005 (2006). [CrossRef]
25. C. G. Schoer, “Reconstructing X-ray fluorescence microtomograms,” Appl. Phys. Lett. 79, 1912 (2001). [CrossRef]
26. B. Golosio, A. Simionovici, A. Somogyi, L. Lemelle, M. Chukalina, and A. Brunetti, “Internal elemental micro-analysis combining x-ray fluorescence, Compton and transmission tomography,” J. Appl. Phys. 94, 145 (2003). [CrossRef]
27. P. J. Green, “On the use of the EM algorithm for penalized likelihood estimation,” J. R. Stat. Soc. 52, 443–452, (1990).
28. E. Levitan and G. T. Herman, “A maximum a posteriori probability expectation maximization algorithm for image reconstruction in emission tomography,” IEEE T. Med. Imaging 6, 185–192, (1987). [CrossRef]
29. J. Moussouris, “Gibbs and Markov random systems with constraints,” J. Stat. Phys. 10, 11–33, (2014). [CrossRef]
30. K. Lange, “Convergence of EM image reconstruction algorithms with Gibbs smoothing,” IEEE T. Med. Imaging 9, 439–446, (1990). [CrossRef]
31. J. H. Chang, J. M. M. Anderson, and J. R. Votaw, “Regularized image reconstruction algorithms for positron emission tomography,” IEEE T. Med. Imaging 23, 1165–1175, (2004). [CrossRef]
32. F. De Carlo, D. Gürsoy, F. Marone, M. Rivers, D. Parkinson, F. Khan, N. Schwarz, D. Vine, S. Vogt, S. C. Gleber, S. Narayanan, M. Newville, A. Lanzirotti, Y. Sun, Y. Hong, and C. Jacobsen, “Scientific data exchange: a schema for HDF5-based storage of raw and analyzed data,” J. Synchrotron Radiat. 21, 1224–1230, (2014). [CrossRef] [PubMed]
33. T. Bicer, “Supporting data-intensive scientific computing on bandwidth and space constrained environments”. PhD Dissertation, (2014).
34. D. Gürsoy, F. De Carlo, X. Xiao, and C. Jacobsen, “TomoPy: a framework for the analysis of synchrotron tomographic data,” J. Synchrotron Radiat. 21, 1185–1193, (2014). [CrossRef]
