1. Introduction

The use of characteristic x-ray emission lines to distinguish between different chemical elements in a specimen goes back to the birth of quantum mechanics [1]. X-ray fluorescence can be stimulated by energy transfer from electron or proton beams, but the best combination of sensitivity and minimum radiation damage is provided by using x-ray absorption [2–4] for this purpose. This is usually done in a scanning microscope mode, where a small x-ray beam spot is raster-scanned across the specimen while x-ray photons are collected by an energy dispersive detector that provides a measure of the energy of each emitted photon [5]. Following early demonstrations [6], x-ray fluorescence microscopy is now commonplace in many laboratories and in particular at a wide range of synchrotron radiation light source facilities worldwide. Because the x-ray beam from synchrotron light sources is usually linearly polarized in the horizontal direction, the energy dispersive detector is usually located at a position 90° in the horizontal from the incident beam so as to be centered on the direction of minimum elastic scattering as shown in Fig. 1 (other energy dispersive detector positions can be used [7], with various relative merits [8]). Because the depth of focus of the x-ray beam is usually large compared to the specimen size, one can treat the incident x-ray beam as a pencil beam of constant diameter and thus translate and rotate the specimen to obtain a set of 2D projections from each x-ray fluorescence line for tomographic reconstruction [9,10], even for elements present at low concentration such as trace elements in biological specimens [11].

 

Fig. 1 Top view of the geometry used in x-ray fluorescence microscopy. The x-ray beam is treated as a pencil beam in the z direction that is raster-scanned across the specimen in 1D in the x direction, and the specimen is then rotated before another image is acquired (successive 2D planes are imaged by motion of the 3D specimen in the y direction, into/out of the plane of this top view). The x-ray transmission signal (absorption) is recorded, and the x-ray fluorescence (XRF) signal is recorded over an angular range of Ω_v by using an energy dispersive detector located at 90° to the beam, in the direction of the elastic scattering minimum for a horizontally polarized x-ray beam. The grid overlay on the specimen shows its discretization with a pixel size of L_v; the set of pixels (in 2D; voxels in 3D) through which the XRF signal might undergo self-absorption in the specimen is indicated in orange.

Download Full Size | PDF

A common approach is to collect the photon counts within predetermined energy windows, or to use per-pixel spectral fitting [12], so as to get immediate elemental concentration maps. These maps are then used in conventional tomographic reconstruction schemes, such as filtered back projection and iterative reconstruction techniques [13, 14]. A more recent approach has been to use a penalized maximum likelihood estimation method on the per-pixel spectra recorded by the energy dispersive detector for improved quantification and elemental separation [15]; we refer to this as full-spectrum analysis. A separate complication involves the correction of fluorescence self-absorption, where characteristic x-ray fluorescence photons emitted from one voxel in a 3D specimen might suffer absorption in specimen material that lies between this voxel and the energy dispersive detector. There have been interesting approaches to correct for self-absorption as will be discussed below, but these approaches have not been combined with full-spectrum analysis. For these reasons, we consider here a combined approach that incorporates both full-spectrum fluorescence analysis, and transmission imaging using absorption, as part of an optimization-based approach to fluorescence tomography analysis by using a complete forward model of the x-ray imaging process, and provide an initial demonstration of our approach on a set of experimental data obtained from the Advanced Photon Source.

2. Fluorescence self-absorption

We begin with a simple illustration of self-absorption. Consider a specimen that consists of a 200 μm diameter borosilicate glass cylinder with a 10 μm diameter tungsten wire off to one side, and a 10 μm gold wire off to another side (Fig. 2). The borosilicate glass was assumed to consist of 81% SiO₂, 13% B₂O₃, 3.5% Na₂O, 2% Al₂O₃, and 0.5% K₂O, with a density of 2.23 g/cm³. If a chromatic x-ray focusing optic like a Fresnel zone plate is used to produce the x-ray pencil beam, a monochromatic x-ray beam should be used for scanning and its energy might be set to 12.1 keV to be well-separated in energy from the Au Lβ₁ fluorescence line at 11.4 keV. As the specimen is rotated, one obtains x-ray transmission (XRT) sinograms based on the attenuation of 12.1 keV x rays in Si, W, and Au as shown at right in Fig. 2. The situation with the x-ray fluorescence (XRF) signal is different; when the x-ray beam is at the right edge of the Si cylinder, the Si Kα₁ x-rays with 1.74 keV photon energy will have to traverse nearly 200 μm of Si before they reach the XRF detector located at left (Fig. 1). Since Table 1 shows that the absorption length of 1.74 keV x-rays is 1.66 μm in glass, Si x-ray fluorescence in particular will be strongly self-absorbed. The fraction of the signal reaching the XRF detector is only exp[−200/1.66] ≃ 5 × 10⁻⁵³ so that essentially none of the Si XRF signal is detected in this case. In fact, only the Si XRF signal from the side nearest to the XRF detector is registered, so that from the Si XRF signal one cannot distinguish between a solid Si cylinder versus one that is hollowed out as shown in the bottom row of Fig. 2. The Au and W XRF signals can better traverse through the Si cylinder to reach the XRF detector, and moreover the 12.1 keV incident beam is also only partly absorbed so by combining all of these signals one can indeed distinguish between a solid and hollow Si cylinder.

 

Fig. 2 Illustration of the x-ray fluorescence self-absorption effect, and how x-ray transmission can be used to recognize and correct for it. We show here a specimen composed of cylinders, or circles in this top view. The largest is of borosilicate glass (composition described in Sec. 2) with 200 μm diameter, followed by tungsten (W) with 10 μm diameter, and gold (Au) with 10 μm diameter. As 1D scans in beamlet positions τ are collected at successive specimen rotation angles θ, one builds up sinograms or (τ, θ) views of elemental x-ray fluorescence (XRF) signals such as the Si XRF signal shown in the middle, as well as 12.1 keV x-ray transmission (XRT) sinograms as shown at right (based on absorption contrast). If there is no self-absorption of the fluorescence signal, one obtains a Si XRF sinogram as shown in the top row, where the incident x-ray beam is partially absorbed in the small W and Au wires as they rotate into positions to intercept the incident beam before it reaches the glass cylinder. However, the 200 μm diameter glass cylinder is large compared to the 1.66 μm absorption length μ⁻¹ of Si Kα₁ x-rays in the glass as shown in Table 1, so that a fraction 1 − exp[−200/1.66] ≃ 1 – 5 × 10⁻⁵³ of the Si XRF signal will be self-absorbed in the rod. As a result, the Si XRF signal will be detected mainly when the incident beam is at the left side of the scan; this leads to the Si XRF sinogram shown in the middle row (the sinogram also shows absorption of the Si XRF signal in the W and Au wires as they rotate through positions where they partly obscure the XRF detector’s view of the Si cylinder). In the bottom row we show the case where the glass cylinder is hollow, with a wall thickness of 30 μm that is nevertheless large compared to the 1.66 μm absorption length of Si XRF photons; in this case the Si XRF sinogram is almost unchanged, but the XRT sinogram is clearly different. By using the combined information of the fluorescence (XRF) and transmission (XRT) sinograms, one can in principle obtain a better reconstructed image of the specimen in the case of strong fluorescence self-absorption.

Download Full Size | PDF

Table 1. X-ray absorption lengths μ⁻¹ for silicon in a borosilicate glass, tungsten and gold at the energies of selected x-ray fluorescence lines and 12.1 keV as incident x-ray energy.

View Table | View all tables in this article

Correction for the self-absorption effect can be achieved from either algorithmic or experimental acquisition perspectives. One possibility to minimize data acquisition time is to use a half rotation/full translation scanning method [16]. From an algorithmic perspective, several methods have been used to correct for the self-absorption effect, including earlier approaches used for radionuclide emission tomography [17–19]. If one can measure the transmission sinograms of the specimen at the energies of all x-ray fluorescence lines, it is possible to correct for self absorption [20]. However, this approach is exceedingly difficult to realize experimentally, since a large number of x-ray fluorescence lines are present in many specimens and one would need to collect a transmission tomography dataset at each of these energies. In the case of uniform absorption and illumination at a single x-ray energy, analytical approaches have been developed [21, 22] and these have been shown [14] to provide a good starting point for the iterative methods we now describe. One approach is to use algebraic (rather than filtered back projection) reconstruction methods to better handle limited rotational sampling, and least-squares fitting to better handle quantum noise [23]; other approaches have used ordered-subsets expectation maximization [24]. In more recent work, optimization approaches have been introduced where the transmission tomography data at a single x-ray energy was used to estimate the absorption at all x-ray fluorescence energies using the fact that (in the absence of x-ray absorption edges) x-ray absorption scales in a power-law relationship with x-ray energy [25–27]. One can also add the Compton scattered signal as another measurement of overall specimen electron density, and use the tabulated absorption coefficients ${\mu }_e^E$ of all elements e at each fluorescence energy E [28]. Other approaches classify the specimen as being composed of a finite number of material phases for the calculation of self-absorption [29]. The optimization approaches in particular serve as inspiration for our approach, which we believe is unique in combining both full-spectrum analysis and transmission imaging along with fluorescence.

We begin by briefly describing the mathematical “forward models” of XRF and XRT. Next, we detail our joint reconstruction approach and the formulation of the objective function and corresponding optimization algorithm. We then discuss choices of scaling parameters in the numerical implementation of the algorithm and present the performance of our joint inversion compared with existing approaches on real datasets.

3. Mathematical model

We start from an earlier approach [30], which we extend considerably here to include a different model of XRF self-absorption effect and the ability to better balance differences in variability of acquired data. We use θ ∈ Θ and $\tau \in \mathcal{T}$ to denote, respectively, the index of the x-ray beam angle and discretized beamlet from a collection of |Θ| angles and $|\mathcal{T}|$ beamlets. The set $\mathcal{V}$ denotes the collection of $|\mathcal{V}|$ spatial voxels used to discretize the reconstructed sample. By $\mathbf{L}=[L_v^{\theta ,\tau }]$, we denote the intersection lengths (in cm) of beamlet (θ, τ) with the voxel $v\in \mathcal{V}$. We use ℇ to denote the collection of |ℇ| possible elements and ${\mu }_e^E$ to denote the mass attenuation coefficient (in cm²g⁻¹) of element e at beam incident energy E. Our goal is to recover $\mathcal{W}=[{\mathcal{W}}_{v,e}]$, the concentration (in g cm⁻³) of element e in voxel v.

3.1. XRT imaging model

A conventional way (see, e.g., [31]) to model the XRT projection ${\tilde{F}}_{\theta ,\tau }^T$ (in counts/sec) of a sample from beamlet (θ, τ) is

To better explore the correlation of XRF and XRT and to link these two data modalities by the common variable $\mathcal{W}$, we note that the coefficients ${\tilde{\mathit{\mu }}}^E$ depend on $\mathcal{W}$ by way of ${\tilde{\mu }}_v^E={\displaystyle \sum _e{\mathcal{W}}_{v,e}{\mu }_e^E}$. Incorporating this fact in Eq. (1), we obtain a new XRT forward model based on $\mathcal{W}$

To obtain a simple proportional relationship, we divide both sides of Eq. (2) by I₀ and then take the logarithm to obtain the XRT forward model used in this work:

We similarly take the logarithm of the raw XRT sinograms used in this paper.

3.2. XRF imaging model

Our discrete model follows an elemental approach, in the sense that we model the XRF energy emitted from a single elemental atom by its corresponding elemental unit spectrum. We first apply spectral blurring to each unit spectrum according to the detector’s response function. Then, justified by the fact that photon counts are additive, the total XRF spectrum detected from the given specimen model comes in after which is estimated as a weighted sum of the unit spectra of the elements being recovered.

First, we model the net XRF intensity I_e,ℓ,s, which corresponds to the characteristic XRF energy E_e emitted from element e, by Sherman’s equation [32] up to first order (i.e., neglecting effects such as Rayleigh and Compton scattering):

In our calculations, the quantity ${\omega }_{e,\ell }\left(1-\frac{1}{r_{e,s}}\right){\mu }_e^E$ is approximated by the XRF cross sections provided from xraylib [33]. Next, we convert the intensity to a spectrum by incorporating the practical experimental environment. Given an energy-dispersive XRF detector with energy channels x_i, $i=1,\dots ,|ℐ|$, we define an indicator function

Then we have the ideal, delta-function peak ${\mathbf{I}}_{e,l,s}^{\mathbf{x}}=I_{e,l,s}{\mathbf{1}}_{E_e}^{\mathbf{x}}$. In practice, because of the detector energy resolution [2], discrete x-ray lines get broadened by a Gaussian distribution with a standard deviation σ. The resulting unit spectrum of element e is thus given by ${\mathbf{M}}_e={\displaystyle \sum _{\ell ,s}{\mathbf{M}}_{e,\ell ,s}}$, where

We then model the total XRF spectrum of a sample with multiple elements by explicitly considering the attenuation of the beam energy and the self-absorption effect of the XRF energy. We represent the attenuation experienced by beamlet (θ, τ) (at incident beam energy E) as it travels toward voxel v by

We let $P_{v,e}^{\theta ,\tau }(\mathcal{W})$ be the attenuation of XRF energy emitted from element e at voxel v by beamlet (θ, τ). To reduce the complexity of the calculation, instead of tracking all the emitted photons isotropically, we consider only the emission from the region Ω_v. This region is the part of the sample discretization that intersects the pyramid determined by the centroid of the voxel v and the XRF detector endpoints; see Fig. 1 for a 2D illustration. In a slight abuse of notation, we let v″ ∈ Ω_v indicate that the centroid of voxel v″ is contained in the region Ω_v. Then the self-absorbed XRF energy is approximated by

4. Optimization-based reconstruction formulations and algorithms

We take $D_{\theta ,\tau }^{\tilde{T}}\in ℝ$ and ${\mathbf{D}}_{\theta ,\tau }^R\in {ℝ}^{|ℐ|}$ to denote the experimental data for XRT and XRF, respectively. We now solve reconstruction problems involving the models $F_{\theta ,\tau }^T(\mathcal{W})$ and ${\mathbf{F}}_{\theta ,\tau }^R(\mathcal{W})$. Given that both these data sources are derived from measured photon counts, we follow a maximum likelihood approach that assumes the measurements are subject to independent Poisson noise [35, 36]. First, we take a logarithm of $D_{\theta ,\tau }^{\tilde{T}}$ and work with $D_{\theta ,\tau }^T=-\ln (D_{\theta ,\tau }^{\tilde{T}}/I_0)$. Maximizing the likelihood (derived in App. A) for our joint inverse problem then can be written as

${\tilde{\varphi }}^R$ and ${\tilde{\varphi }}^T$ correspond to the XRF and XRT objective terms, respectively; and β₁ ≥ 0, β₂ ≥ 0 are scaling parameters. The scaling parameter β₁ balances the ability of each modality to fit the data, and β₂ detects the relative variability between the data sources ${\mathbf{D}}_{\theta ,\tau }^R$ and $D_{\theta ,\tau }^T$.

Algorithm 1. Algorithm for Solving Joint Inversion with Linearized XRF Term.

View Table | View all tables in this article

Advances in x-ray sources, optics, and detectors mean that the datasets to be analyzed can be large; thus, having a fast and memory-efficient algorithm to solve (8) is highly desirable. Therefore, we apply an alternating direction approach described in Alg. 1. In this approach, instead of directly minimizing Eq. (8), we first solve an “inner iteration” subproblem:

To approximately solve Eq. (9), we adapt a form of the inexact truncated Newton (TN) method in [37]. We write TN as a function of the form

Algorithm 2. Backtracking Line Search

View Table | View all tables in this article

5. Experimental reconstruction

We now demonstrate the benefit of the proposed joint inversion approach by using experimental data. We constructed a simple test sample consisting of a borosilica glass rod with a composition as described in Sec. 2, wrapped with a W wire of 10 μm diameter, and a Au wire of 10 μm diameter. This test specimen was scanned by an incident beam energy of 12.1 keV at beamline 2-ID-E at the Advanced Photon Source at Argonne National Laboratory. Each projection was acquired by raster scanning horizontally and vertically with a 200 nm scanning step size and 73 scanning angles over an angular range of 360°, with the data stored as HDF5 files with metadata included [38]. At each scanning step, the transmission signal was acquired using a four-segment charge-integrating silicon drift detector [39], while the fluorescence signals were collected by using an energy-dispersive detector (Vortex ME-4) located at 90° relative to the incident beam, covering an energy range of 0–20 keV with 2,000 energy channels. We then sum the signal from the 4 transmission detector elements without distinguishing the fact that each element sees a slightly different angle to obtain the final data; this results in 73 projections, with each slice involving 1, 750 × 51 pixels, and leads to a dataset of dimension 73 × 1, 750 × 51 × 2, 000.

As expected from the attenuation lengths given in Table 1 and the simulations shown in Fig. 2, this dataset shows strong self-absorption in the Si fluorescence measurements. We compare our reconstruction result with the output of TomoPy 0.1.15 [40], a widely used tomographic data-processing and image reconstruction library. TomoPy takes the elemental concentration map decomposed from the raw spectrum by the program MAPS 1.2 [12] for improved photon statistics compared with the raw data. The MAPS program fits the full energy spectrum recorded at each scan to a set of x-ray fluorescence peaks plus background signals, and it returns a 2D dataset corresponding to a certain elemental concentration per scan. Figure 3 shows three elemental XRF sinograms of interest as calculated by this approach. Within TomoPy, we used a maximum likelihood expectation maximization algorithm to reconstruct the three sinograms separately. All numerical experiments are performed on a platform with 1.5 TB DDR3 memory and four Intel E7-4820 Xeon CPUs.

 

Fig. 3 Relative elemental concentrations obtained from a MAPS-based fit of the raw x-ray fluorescence data for the glass rod sample. Due to the imperfection of fitting and background rejection (which might be able to be corrected with additional expert input), the decomposed elemental concentrations show certain artifacts, where certain elemental sinograms pick up other elements’ signals. For example, according to our knowledge of the sample, we know that Si exists only in the rod part with a cylinder shape; but its corresponding sinogram shows that it also exists in the two wires, which is caused by imperfect data fitting. Those two extra curves in the sinogram are actually picked up from Au and W signals because certain emission lines of Au and W overlap those of Si. Similar artifacts happen to Au and W sinograms as well.

Download Full Size | PDF

For the purposes of algorithm testing with reduced computational cost, we reconstructed only a 2D middle slice (x, y) of the 3D glass rod dataset (x, y, z). For each angle, we summed together 9 adjacent y slices of both XRF spectra and XRT measurements as the input experiment data. Furthermore, instead of using the full 2000 energy channels, we pick 40 energy channels around each of the 20 emission lines provided by xraylib for each element’s K_α, K_β, L_α, L_β, and M_α lines; this process returns a total of 800 energy channels to be considered in the reconstruction. Therefore, in our illustration, |θ| = 73, |τ| = 195, |ℇ| = 800, and $|\mathcal{V}|=(|\tau |)\times (|\tau |)$.

5.1. Selection of β values

Next, we explain one way to select values for β₁ and β₂ for use in the objective in Eq. (8). We recall that the effect of β₂ is to balance the magnitudes of the XRF and XRT measurement data, and that its exact value is not critical. Therefore, according to the magnitude difference of two data sources shown in Fig. 4, we chose to use β₂ = 100 to balance this difference so that both measurements have maxima near 15 in their respective units. The selection of β₁ is accomplished by applying the so-called “L-curve method” [41]. In Fig. 5, we plot the L-curve defined by the curve of XRT terms ${\tilde{\varphi }}^T$ versus XRF terms ${\tilde{\varphi }}^R$ obtained from solving Eq. (8) with different β₁ values. This curve displays the tradeoffs between these two modalities and provides an aid in choosing an appropriate balancing parameter β₁. The curvature, defined as the curvature of a circle drawn through three successive points on the L-curve, is calculated and plotted in Fig. 6. As suggested by [41], we choose the point on the L-curve with the maximum curvature; according to the objective values of ${\tilde{\varphi }}^R$ and ${\tilde{\varphi }}^T$ shown in Fig. 5 and the curvature shown in Fig. 6, this is β₁ = 1. Figure 10 shows the reconstructed elemental maps corresponding to different β₁ values. For the two extreme cases β₁ = 0.001 and β₁ = 100, the reconstruction results resemble the single-modality reconstructions. However, there is no significant difference between the results returned for β₁ ∈ [0.25, 10] based on our prior knowledge of the sample. This lack of sensitivity to precise balancing parameter values indicates robustness for practical utilization.

 

Fig. 4 Experimental sinograms. Left: mean (across energy channels) value of XRF raw spectrum; Right: XRT optical density. Based on the different magnitudes of these two datasets, we choose β₂ = 100 as the scaling parameter to balance the measurement variability of the two data sources, so that both measurements have maxima near 15 in their respective units. As a result, the relative variability of the two detectors between the data sources is mitigated.

Download Full Size | PDF

 

Fig. 5 Method for choosing the parameter β₁ that appears in the cost function of Eq. (8): XRF objective value ${\tilde{\varphi }}^R$ versus XRT objective value ${\tilde{\varphi }}^T$ given different values β₁, with fixed β₂ = 100. The curve displays the tradeoff between these two modalities.

Download Full Size | PDF

 

Fig. 6 Method for choosing the parameter β₁ that appears in the cost function of Eq. (8): The curvature from successive points in Fig. 5; the point with maximum curvature occurs at β₁ = 1.

Download Full Size | PDF

5.2. Joint inversion results

Given ${\mathcal{W}}^0=\mathbf{0}$ as the initial guess for the joint inversion, Fig. 7 shows the reconstruction result for each element by using Alg. 1 with ϵ = 10⁻⁶. In particular, Fig. 8 shows the performance of the inner iteration by TN to reduce both the XRF and XRT objective values. Correspondingly, Fig. 9 shows the reconstructed result of each outer iteration of Alg. 1. The reconstructed elemental maps show the benefits of our joint inversion mainly from two perspectives. First, because of the imperfections of spectral fitting and background rejection, the decomposed elemental concentrations show certain artifacts—which we call the “elemental crosstalk ”—where certain elemental sinograms pick up other elements’ signals. For example, according to our knowledge of the sample, we know that Si exists only in the rod part with a cylinder shape; but its corresponding sinogram from MAPS (Fig. 3) would suggest that it also exists in the two wires; this is caused by imperfect data fitting. Those two extra curves in the sinogram are actually picked up from Au and W signals because certain emission lines of Au and W overlap those of Si. As a result, the reconstruction from TomoPy based on these decomposed sinograms will contain the “elemental-crosstalk points,” which are shown as two small dots around Si and a dot in the W map in the left bottom corner of Fig. 7. Comparing the results from XRF single inversion using our forward model with the TomoPy output, we see that our forward model is able to better distinguish the different elemental signals; that is, the “elemental crosstalk” is greatly mitigated. Furthermore, by introducing the XRT modality into the reconstruction, the joint inversion not only suppressed the artifacts from the “elemental crosstalk” introduced by the imperfect fitting and background rejection, but also more accurately recovered Si by filling the inside of the cylinder and thereby correcting the self-absorption effect. Also, we provide a quantity evaluation of our reconstruction. We simulate the XRF spectrum based on the reconstructed elemental composition and compare it to the real experimental data as shown in Fig. 11. We can see that, except the background region that our forward model does not include to simulate, the essential peaks corresponding to the main elements we are interested to cover agree very well with the experimental data. This comparison indicates that not only our joint reconstruction improves the solution from a visualization perspective, but also from quantification point of view. Furthermore, it indicates a satisfying accuracy of our XRF forward model.

 

Fig. 7 Comparison of reconstruction results for MAPS+TomoPy, XRF alone, and joint reconstruction, respectively, given an initial guess of all zeros, β₁ = 1, and β₂ = 100. Every elemental map is rescaled to the range [0, 0.5]. It is clear that the joint reconstruction returns the best result from two perspectives: first, the glass rod is filled with Si; and second, the “elemental crosstalk” is dramatically mitigated for the reconstruction of Si and W.

Download Full Size | PDF

 

Fig. 8 Convergence of TN for each inner iteration j, given a maximum number of inner iterations as 52, β₁ = 1, and β₂ = 100. We can see that along the iterations, TN is reducing the objective function so that the forward model fits better and better to the given data.

Download Full Size | PDF

 

Fig. 9 Solution for each outer iteration i, given (β₁, β₂) = (1, 100). At iteration i = 3, Alg. 1 reaches its stopping criterion in the sense that the solution does not change anymore. The results also indicate that our alternating algorithm requires very few outer iterations.

Download Full Size | PDF

 

Fig. 10 Reconstruction results (given β₂ = 100) for different β₁ values. Apart from the two extreme cases (β₁ = 0.001, where XRF dominates, and β₁ = 100, where XRT dominates), the reconstructions do not show clear difference in terms of quality. Therefore, for a broad range of β₁ values, the joint reconstruction is able to dramatically improve upon the single-modality reconstructions.

Download Full Size | PDF

 

Fig. 11 Example x-ray fluorescence spectrum. In this case, an incident beam with a photon energy of 12.1 keV is used to excite x-ray fluorescence from a specimen consisting of a borosilicate glass cylinder comprised mainly of SiO₂ but with other elements present, and tungsten (W) and gold (Au) wires. The experimental spectrum is averaged over all positions of a sinogram from one scan row. The simulated spectrum based on the reconstructed elemental map is generated by the forward model described in Sec. 3; it includes tabulated [33] x-ray fluorescence lines for all elements present in the specimen along with the Gaussian energy response of the fluorescence detector, plus the background spectrum from non-specimen areas. Some additional background is present in the 4–7 keV energy range due to the materials in the experimental apparatus as indicated at specific fluorescence peaks; because this background does not change whether or not a specimen region is illuminated, it does not affect our analysis.

Download Full Size | PDF

One limitation of our method is the high requirement of computing resources for generating the forward mapping tensor incorporating the known mass attenuation coefficients. Overall, the reconstruction time spent on this particular test is 4301 seconds plus a couple of hours on generating the forward mapping. The main bottleneck is in the large memory requirement and the calculation of the self-absorption term. Parallelizing our TN solver and projections/beamlets, or using multilevel schemes [42], would also accelerate the performance. This will be the focus of our future work.

6. Conclusion

Guided by the multimodal analysis methodology developed in [30], we apply a joint-inversion framework to solve XRF reconstruction problem more accurately by incorporating a second data modality as XRT. We investigate the correlations between XRF and XRT data, and establish a link between datasets by reformulating their models so that they share a common set of unknown variables. We develop an iterative algorithm by alternatively maximizing a Poisson likelihood objective to estimate the unknown elemental distribution, and then updating the self-absorption term in the forward model. The initial demonstration presented in the paper show that when facing strong self-absorption effects, significant improvements are achieved by performing joint inversion. Furthermore, because of the improved accuracy provided by our XRF forward model, the artifacts arising from the “elemental crosstalk” are greatly mitigated.

The bottleneck of the current code is in its extensive memory requirement to evaluate large-scale tensor product involving raw spectra with many energy channels. In the future, we expect that larger-size problems will be achievable once we move beyond our prototype code. Parallelizing our TN solver and projections/beamlines, or using multilevel schemes [42] would also accelerate the performance. On the other hand, we will explore the combination of our method with different type of data acquisition (e.g., suggested in [16]) to achieve better performance.