Quantitative X-ray fluorescence computed tomography for low-Z samples using an iterative absorption correction algorithm

X-ray fluorescence computed tomography is often used to measure trace element distributions within low-Z samples, using algorithms capable of X-ray absorption correction when sample self-absorption is not negligible. Its reconstruction is more complicated compared to transmission tomography, and therefore not widely used. We describe in this paper a very practical iterative method that uses widely available transmission tomography reconstruction software for fluorescence tomography. With this method, sample self-absorption can be corrected not only for the absorption within the measured layer but also for the absorption by material beyond that layer. By combining tomography with analysis for scanning X-ray fluorescence microscopy, absolute concentrations of trace elements can be obtained. By using widely shared software, we not only minimized the coding, took advantage of computing efficiency of fast Fourier transform in transmission tomography software, but also thereby accessed well-developed data proce...


I. INTRODUCTION
X-Ray fluorescence computed tomography (XFCT) has been used in the synchrotron community for decades, 1-7 commonly implemented with a combination of translation and rotation scans using pencil beam X-rays.If X-ray absorption by the sample is negligible, XFCT reconstruction can be performed in the same way as transmission X-ray computed tomography (CT), [8][9][10][11][12] for which many well-developed software packages are available, either as open source or commercial product.However self-absorption is generally not negligible and many reconstruction methods have been developed to correct such absorptions.These methods include algebraic reconstruction techniques (ART), [13][14][15][16] ART with ordering scheme, 17 maximum likelihood expectation maximization (ML-EM), 15,18 ordered subset EM algorithm (OS-EM), 19,20 least-squares method using singular value decomposition, 21 reconstruction with an approximation proposed by Hogan et al, [22][23][24] penalizedlikelihood reconstruction, 25 conjugate gradient method, 26 and iterative refinement. 27Compared to transmission CT reconstruction programmed with fast Fourier transform (FFT), other algorithms with self-absorption correction are more complicated and considerably more time consuming, and therefore not widely used.
Meanwhile, to simplify the work, there have been attempts to quantify XFCT using an absorption correction on the sinograms (the 2D array of projection data) without considering the distributions of fluorescent elements in the sample. 28,29Strictly speaking, such absorption correction is only an approximation, or a correction to the first order.Experimentally, a setup with dual X-ray fluorescence (XRF) detectors has been used to mitigate absorption to certain degree. 302][33][34][35] Depending on both the sample sizes and the XRF attenuation lengths of measured elements, it is often necessary to correct the self-absorption during XFCT reconstruction.The main composition of such samples is low Z material with fluorescence energy too low to excite secondary fluorescence of the measured trace elements.For these samples it is possible to have the self-absorption corrected in a simple way as presented in this paper.We will demonstrate that, for trace element analyses in such low-Z samples, XFCT reconstruction can be performed using transmission CT software with iterative absorption correction on sinograms.That makes XFCT reconstruction easier because tools for transmission CT are widely available, much faster, well developed and supported.In general, the absorption correction of existing algorithms is limited to absorption in the measured layer.With this new algorithm, for each measured layer, the self-absorption by a sample's three-dimensional structure can also be corrected.At CHESS, a Python user interface called Praxes has been developed using PyMca libraries 36 for scanning X-ray fluorescence microscopy (SXFM) data processing.For quantitative XFCT for biological samples, absolute concentration can be calculated by data pre-processing with Praxes before the reconstruction.By using widely shared and well-developed software packages, we not only minimized the coding work for self-absorption correction, but also gained access to all the widely used data processing tools coming with these software packages for reliable, high quality work.Therefore we suggest that our reconstruction approach is a very practical way for XFCT with absorption correction for low-Z samples.
Fish otoliths (literally, "ear stones" that form part of the hearing and balance system in fishes) and eye lenses grow incrementally throughout the lifetime of fishes, and thus have potential as archival structures.8][39] For fish otoliths, the calcium carbonate material is hard enough to be thin sectioned; its high density also limits X-ray penetration and helps preserve spatial resolution with our ordinary synchrotron SXFM setup, which orients the sample surface 45°to both the incident beam and the fluorescence detector.In contrast, slicing a proteinaceous fish eye lens to about 100 µm thin or less is not always successful, as the softer material tends to flake off or break.When the lens section is thick enough not to break apart, high penetration of X-rays through the low-Z material of lens compromises the spatial resolution of SXFM 2D mapping (i.e., elements fluoresce below the surface and are imaged).Therefore, XFCT of an unsliced eye lens for trace element measurement seems very promising, if XFCT self-absorption can be corrected.As an example, the XFCT reconstruction of a fish eye lens's absolute concentrations of several trace elements is demonstrated in this paper.The results are quantitatively compared with the XRF 2D mapping of the sliced lens of the same fish's other eye, showing not only a good agreement, but also a better spatial resolution with XFCT.

A. XFCT data pre-processing
The experimental setup of XFCT is similar to SXFM, with an addition of sample angular rotation.The open source Python software package Praxes (https://github.com/praxes/praxes),developed at CHESS by Dr. Darren Dale for SXFM, is used for XFCT data pre-processing.Praxes is a package that wraps PyMca libraries 36 to handle XRF spectrum fitting, calibration and calculation of element absolute concentration.It can run in real-time mode, where the execution of XRF scans and real-time data analysis can take place simultaneously.It can also be used for data processing after an XRF scan is completed.Currently, this software utilizes the Cython module (C extensions for Python) and the Python multiprocessing module to speed up calculation-intensive, time consuming processes.
Trace element concentrations are calculated by PyMca libraries using the fundamental parameter method and the results are presented in units of mass fractions.PyMca assumes a layered structure of samples.Within each pixel, the trace element distribution is assumed to be uniform along the X-ray path in each layer.Since such samples are mainly made of low Z elements and their fluorescence energy is too low to excite trace elements, secondary fluorescence is negligible and ignored in calling PyMca libraries.Currently only the fundamental parameter method in PyMca is implemented through Praxes.In pre-processing we treat the three-dimensional sample as a very thin sheet of layered structure, therefore X-ray absorption is not considered.The correction of self-absorption is only handled during tomographic reconstruction.The calibration of experimental setup is performed using either an otolith standard or a pure and thin metal foil.The XFCT data after pre-processing with Praxes provides sinograms in units of element concentration integrated along the incident X-ray beam path in the sample.A Vortex ME4 silicon drift detector (SDD) is used for XFCT measurement.

B. Iterative absorption correction on XFCT sinograms
The difference between X-ray transmission CT and XFCT can be explained with the following equations: where I(θ,t) is the logarithm of the measured sinogram for transmission CT, and I f (θ,t) the sinogram for XFCT; p(x,y) is the X-ray linear attenuation coefficient for transmission CT, and p(θ,s,t) is proportional to element concentration for XFCT.While the (x,y) axes are fixed in the reference frame of the sample, the (s,t) axes are always parallel to the X-ray beam and sample translation at any rotation angle θ, therefore, s = xcosθ + ysinθ, t = xsinθ + ycosθ.Eq.( 1) is known as the Radon transform and represents transmission CT projections. 14Many reconstruction software packages are available to recover p(x,y) from measured I(θ,t).For XFCT as described by eq.(2), f(θ,s,t) represents the transmission of incident X-rays along the beam path up to position (s,t) inside the sample, and g(θ,s,t) represents the integrated attenuation of XRF by the sample from position (s,t) up to the XRF detector. 15,20,23,26The values of f, g at each voxel (a grid in 3D space) (θ, s, t) are related to element distributions within the sample.If we re-write eq.(2) as: or, and if we call T(θ, t) the absorption factor, at each point (θ, t), this absorption factor depends not only on the fluorescent element distribution and sample major composition distribution along the incident beam path, but also on the sample composition distribution along all the XRF paths to the detector.That makes the direct reconstruction of XFCT complicated.On the contrary, if the distributions of sample composition are known, it is straightforward to calculate the absorption factor T(θ,t); with measured I f (θ, t) and calculated T(θ,t), the XFCT reconstruction of eq.( 4) becomes equivalent to the transmission CT of eq.(1).For many biological samples, the major composition is approximately known, and the sample's 3D structure can be reconstructed using the transmitted X-rays measured simultaneously during XFCT scans.However, the fluorescent trace element distribution needed for the absorption factor calculation is initially unknown and is actually what we would like to measure.Therefore it seems appropriate to utilize an iterative method: (i) to calculate absorption factor T with assumed trace element distribution p(x,y); (ii) to reconstruct p(x,y) according to eq.( 4) using measured I f (θ,t) and calculated T(θ,t); (iii) to repeat (i) and (ii) until a converged solution is achieved.X-ray attenuation coefficients at both incident energy (µ 0 ) and fluorescent energy (µ i ) are needed to calculate the absorption factor.While µ 0 can be directly measured using transmitted X-rays during XFCT scans, µ i needs to be calculated either according to sample composition 17,40 or by utilizing the sinogram asymmetric information arising from self-absorption. 26In this paper we use sinogram asymmetry for the µ i calculation and absorption corrections.We use µ i calculated from the sample composition as a reference for cross-checking.Similar work was done using the conjugate gradient method under physical considerations for a µ i calculation and conjugate gradient method for reconstruction. 26owever, we calculate an average µ i in each iteration just simply using the slope of the asymmetric sinograms d(log(I fc (θ,t)/I fc (θ+π,t)))/dt, and using transmission CT software for the reconstruction.In this way it is not only possible to correct self-absorption by material beyond the measured layer using sample 3-dimensional structure, but also it minimizes the coding work by using widely available transmission CT software.
To summarize, the complete procedure of quantitative XFCT reconstruction is: (1) use transmission CT reconstruction (such as filtered back-projection) to obtain the sample structure, and use the same CT reconstruction method to obtain an initial element distribution p(x,y) after the XFCT data are pre-processed with Praxes; (2) initialize µ i according to known (or assumed) sample major composition; (3) calculate the absorption factor T(θ,t) using the previously calculated element distribution, sample 3D structure, and information of 4-element Vortex detector positions.(4) Use the filtered back-projection to obtain p(x,y) according to eq.( 4), followed by updating T(θ,t) (or I fc ) and µ i with new p(x,y), and then update T(θ,t) (or I fc ) again with the updated µ i .Repeat steps (3) and (4) until a converged solution is achieved.Here, the open source software TomoPy, 41 including its pre-reconstruction and post-reconstruction data processing tools, is used for CT reconstruction.By using well-established SXFM software and CT software, additional coding is only needed for the straightforward absorption factor calculation, and some script to link the iterations, therefore it is easy to implement.The calibration of the experimental setup for quantitative purposes is identical to that for SXFM.

III. XFCT RESULTS OF FISH EYE LENSES
The diameters of fish eye lenses can be one to several millimeters.The lens is made of crystalline proteins.It is assumed to have a composition similar to the nominal data of the NIST eye lens standard (ICRU-44) (http://physics.nist.gov/PhysRefData/XrayMassCoef/tab2.html), of which the mass fraction by weight is, H: 0.096, C: 0.195, N: 0.057, O: 0.646, S: 0.003, and with density of 1.07 g/cm 3 .The sensor area of the Vortex ME4 SDD is about 180 mm 2 and the separation between the 4 elements is about 12 mm.Depending on the size of the eye lens and the fluorescence intensity level, the distance from the SDD to the sample can be as short as 30 mm, which makes it necessary to correct for self-absorption using the spherical-shaped lens 3D structure and to take into account the geometry of SDD sensor position and size.The experiment was performed at the CHESS F3 bending magnet beamline with a 0.6% bandwidth multilayer monochromator and with a 20 µm beam focused by a glass capillary, 42 with a typical flux above 5x10 10 ph/s at energy around 16-17 keV.

A. Simulation of absorption correction
Simulation over a spherical phantom of a 4.0 mm diameter fish eye lens shows that this iterative method does work.To better represent the lens CT experiment, multiple ellipses were used to mimic the circular distribution of a trace element on a lens cross-section (Fig. 1a).XFCT sinograms were simulated with consideration of absorptions of both incident X-rays (assumed as 16.5 keV) and fluorescence.With the above mentioned lens composition and density, Fe, Zn and Br fluorescence attenuation lengths in the simulated lens (defined as the penetration depth where X-ray intensity decreases to 1/e of its initial intensity) were 0.54 mm, 1.36 mm and 3.65 mm respectively, and the 16.5 keV incident beam attenuation length was 9.3 mm.The reconstructions of these elements without absorption correction (Fig. 1b-1d) would underestimate the intensity of element concentration, especially towards the lens center and for lower Z elements.With iterative absorption correction on the sinograms, the original pattern was recovered (Fig. 1(f)), only showing obvious non perfection in the extreme case of Fe where X-ray attenuation was too severe (Fig. 1(e)).Given the Fe fluorescence attenuation length of 0.54 mm, and sample diameter of 4 mm, Fe fluorescence from the sample center was severely attenuated in CT projections, causing poor reconstruction.However even in this extreme case, the absorption correction provided tremendous improvement on results compared to non-corrected reconstruction (Fig. 1(b)).Simulations of absorption corrected XFCT for Ni, Cu, Se, Br all have generated good results that are visually similar to Fig. 1(a) and Fig. 1(f).To test the capability of µ i search of this method, the initial attenuation length was intentionally set quite far from its true value, either 0.5 or 1.5 times the actual attenuation length of each element.The accurate attenuation coefficient of incident X-ray µ 0 was used in all simulations because in real XFCT experiments it can be measured simultaneously with transmission CT.Fig. 2a shows that it only takes a couple of iterations for µ i to approach a stable value, similar to the observation by a previously published method. 26For fluorescence that is not severely attenuated, the reconstruction converges quickly as µ i converges.However, for the difficult cases where fluorescence is severely attenuated such as exemplified by Fe above, more iterations are needed for the reconstruction to converge (Fig. 2b).To ensure accuracy, we used other converging criteria in addition to µ i .One is the root mean square (RMS) variation of corrected sinograms: where q is the iteration number, and N is the total projection number over all the angles and translations.By the end of our simulation, as shown in Fig. 2b, P q was about 2.1×10 -5 for Br, and about 6.4×10 -3 for Fe.Another criterion is the change between iterations of variance of reconstructed results: where V(f q (x,y)) is the variance of reconstructed results in q th iteration (V = (f − f ) 2 ).A minimum or a stable value of δ q around 1x10 -3 to 1x10 -2 was achieved before the iterations ended, also an indication of a converged calculation as explained by previous studies. 43,44The iteration numbers for the simulations are between 4 to 10.

B. XFCT of fish eye lens and the comparison with its 2D scans
This XFCT reconstruction method was used for trace element mapping of the eye lens of a European fish, Common Bream (Abramis brama).The fish was sacrificed, and the lenses (diameters just under 2 mm) were extracted and allowed to air dry.One of the lens pair was subsequently sliced to about 280 um thick, and quantified with 2D XRF mapping, while its mate was used for CT scan along its middle plane.In both cases, a 20 µm FWHM incident X-ray beam focused by a capillary was used; the incident X-ray energy was 16.2 keV.An in-house standard, made by crushing otoliths into a fine powder and pressing into a pellet, and a 15.9 ± 5% µg/cm 2 Cu metal foil were used for the calibration of experiment setup.For 2D XRF mapping, the sliced lens sample was mounted 45°t owards both incident X-rays and the SDD which was 90°to the direct beam for minimal scattered signals while using polarized synchrotron X-rays.For this specific 2D mapping, step size was 18 µm.For XFCT, a 20 µm translation step size and 1°angular step size were used, and with full 360°rotation.As mentioned before, the absolute concentration of 2D XRF mapping was obtained with Praxes.The sinograms of XFCT, in units of concentration integrated over the X-ray path, were also pre-processed with Praxes.Data from an ion-chamber downstream of the sample were used for transmission CT.Reconstruction of transmission CT provided sample "density" information, as shown in Fig. 3, in units of "normalized density" -the ratio of the X-ray attention coefficient obtained from CT over the theoretical attenuation coefficient calculated for the NIST eye lens standard (ICRU-44).Fig. 3 shows that the actual X-ray attenuation coefficient (or "density") is about 50% more than the "standard", but with little variation over different positions.Using the measured "density" in the absorption correction can help correct assumption errors about either sample density or sample composition.The 3D distribution of the density is needed even just for CT reconstruction of one layer; however for this specific case, since the eye lens was more or less spherical, for simplicity, we assumed lens cross-section shape at z distance away from middle plane was scaled down from the middle plane according to r 2 0 − z 2 /r 0 , where r 0 was the average radius on lens middle plane measured from transmission CT.
The theoretical attenuation coefficient of each fluorescence, based on the eye lens composition of the NIST standard, was used as the initial value of µ i for iterative absorption correction.The µ i reached 055111-7 Huang, Limburg, and Rohtla AIP Advances 7, 055111 (2017) FIG. 3. Distribution of X-ray attenuation coefficient on lens middle plane obtained from transmission CT, normalized to the theoretic attenuation coefficient of standard (referred as "normalized density" or "density" in paper).The discrepancy from 1 could come from inaccurately assumed material composition and/or density.The irregularity of the shape is due to the way the lens dried.In the live fish, the lens is spherical.
a stable value very quickly (Fig 4) in the iteration.The RMS variation of the corrected sinograms between different loops dropped quickly (Fig. 5a), and the change of variance of reconstructed results also approached a minimum or stable value around 0.01 or less (Fig. 5b), an indication of convergence according to previous studies. 43,44Only 3-5 iterations were used to reach a stable solution, depending mostly on how severe the fluorescence attenuation was.Sinograms before and after the absorption correction of Fe, Zn and Br (Fig. 6) clearly showed improved symmetry after the correction.Selected this boundary is quite indistinct in the XRF 2D mapping (compare Figs. 7 and 8).The boundary of an inner sphere also appears clearly on Zn, Br and Rb XFCT, but is only vaguely shown in the corresponding 2D mappings.Ni and Zn XFCT also show high concentration just along the surface of eye lens, but that is not as clearly shown in XRF 2D mapping.The relatively poor spatial resolution on 2D mapping comes from the combination of a thick sample (280 µm in this case) and XRF signal averaging along beam path in sample.Another advantage of XFCT is that it is less affected by contamination.Comparing Fig. 7(a) with Fig. 8(a), while Fe is only detected at some spots on lens surface by XFCT, it is observed on random positions of 2D mapping, which may partly come from handling contamination.
The converged values of the attenuation coefficient µ i in Fig. 4 are off a bit from the theoretical values calculated with nominal compositions of the NIST eye lens standard.Among the 6 elements studied, for half of them the µ i converged to a value only a few percent different from what expected for the standard composition, and all but one converged to values less than 20% different from what expected.That error could mainly come from the composition difference between fish eye lens and NIST eye lens standard, recalling that the measured X-ray attenuation was about 50% higher than expected.However the largest deviation of µ i shown in Fig. 4 is for the element Fe, of which the fluorescence signal is overall very weak except at some isolated spots on the lens surface.As shown in Fig. 6, the sinogram of Fe fluorescence was made of mainly a few fine lines on top of an otherwise very low signal.The reconstruction of this weak Fe sinogram signal may be affected by measurement background which could also result in a lower fitted µ i .In short, we do not exclude the possibility of background influence on the reconstruction of very weak XFCT signals.

IV. CONCLUSION
A simple and very practical XFCT reconstruction algorithm using widely-available transmission CT tools has been introduced.Self-absorption on sinograms can be corrected using iterations of a straightforward calculation of X-ray absorption factors.X-ray fluorescence absorption in a sample beyond the measured layer can also be corrected with this method.Combined with SXFM data processing, this algorithm provides a straightforward method to obtain absolute concentration from XFCT of biological samples.With this method, very little coding is required for researchers having access to transmission CT software and SXFM analysis software.Using transmission CT software will let us benefit from rapid reconstruction by implementing FFT in its algorithm.Many well-developed data processing tools coming with these well-known and reliable software packages will help us for high quality work without our "reinventing the wheel" in coding.As an example, the absolute concentration of the trace elements in an eye lens of a fish has been performed with this method.The reconstruction converges faster for less severely attenuated fluorescence.It takes only 3 -5 iterations to reconstruct the 6 selected elements of the specific lens sample, while it takes 4 -10 iterations for the simulations of a larger eye lens with more absorption.Absolute XFCT results from fish eye lens are compared with 2D XRF scans on the sliced sample, with reasonable quantitative agreement.However there is a clear improvement of spatial resolution by using XFCT, because XFCT measurement is not compromised by sample thickness, a problem we have encountered with 2D mapping.In principle this method is not limited to any specific reconstruction algorithm or specific software of transmission CT, therefore it is possible for us to select a CT algorithm that best fits the specific experiment.
For example, we can select reconstruction tools suitable for incomplete projections for samples with incomplete angular access, etc.Although we have only focused our discussions on biological samples so far, this method should also be suitable for all such applications where the absorption is dominated by relatively lower Z elements, such as the correlative imaging study of catalyst metal poisoning. 29

FIG. 1
FIG. 1.(a) Phantom with eccentric ellipses mimicking the ring patterns of trace element distributions inside a 4 mm diameter fish eye lens.(b),(c),(d) are the reconstructions without absorption correction using simulated sinograms of Fe, Zn and Br.(e), (f) are the reconstruction of Fe, Zn patterns with absorption corrected.Obvious imperfection is only noticeable for the extreme case of Fe of which the fluorescence from sample center can be hardly detected because of absorption.

FIG. 2 .
FIG. 2. (a) Converging of attenuation lengths relative to their theoretic values (shown on legend in millimeter), for simulated elements.The initial lengths are set intentionally either 0.5 times or 1.5 times of theoretic numbers to test the converging capability of simulations.(b) RMS variation of corrected sinograms between iterations, showing faster converging for less attenuated fluorescence and slower converging for severely attenuated fluorescence.

FIG. 4 .FIG. 6 .
FIG. 4. Converging of attenuation coefficients versus iterations, for different fluorescences, normalized by their theoretic values for standard.All simulations start with the theoretic attenuation coefficients of the fluorescence in standard.