Nondestructive Quantification of Heavy Elements Through the Analysis of Beam Hardening Artifacts Using Double-Exposure X-ray Computed Tomography: A Theoretical Consideration

We present a theoretical background for heavy element quantification through the intensive analysis of beam hardening (cupping artifacts) in X-ray computed tomography (CT) images. Cupping artifacts resulting from X-ray CT using a polychromatic X-ray source are quantitatively analyzed with an analytical solution for a cylindrical sample of a homogeneous aqueous solution/suspension containing a heavy element. The theoretical solution reveals that the severity of cupping artifacts is strongly dependent on the sample chemistry and the acceleration voltage of the X-ray tube. Careful analysis of this dependency enabled simultaneous determination of the atomic number and molar concentration of the heavy element within a particular estimation error range. Significant improvement in terms of the accuracy of determining the atomic number was achieved by employing double-exposure X-ray CT.


Introduction
Quantification of heavy elements in samples is one of the most important subjects in analytical chemistry [1]. Substances that contain heavy elements are often toxic and hazardous; therefore, nondestructive methods to safely analyze such heavy elements confined in containers are desirable, such as in cargo security [2] and in analysis of polluted wastes, water and soils [3][4][5][6]. The X-ray computed tomography (CT) is one of useful methods for the nondestructive discrimination of heavy elements in samples [7]. However, heavy elements yield undesirable beam hardening artifacts (cupping artifacts) in the reconstructed CT images due to the polychromatic X-ray source [8]. This cupping artifact occurs as a characteristic cupped shape for the profile of the voxel values across a homogeneous sample. The cupping artifact yielding the inhomogeneous profiles even for homogeneous samples could reduce the accuracy of heavy element discrimination when the conventional X-ray CT method is employed. Nakashima and Nakano [9,10] proposed a novel image analysis method using single-exposure X-ray CT (i.e., using a single acceleration voltage for the X-ray tube). The method enables simultaneous estimation of the atomic number and concentration of a heavy element confined in a sample by carefully analyzing the cupping profiles. This method [9,10] can be implemented using conventional CT systems without energy-resolved [11,12] special detectors.
However, the constraints used for the material differentiation provided by the single-exposure CT method were not complete. As a result, this method failed to eliminate heavyelement candidates with abnormally high atomic numbers (e.g., 75 Re and 76 Os) relative to the ground truth (e.g., 58 Ce) [9]. In order to eliminate such candidates, we propose the use of double-exposure X-ray CT (i.e., acquisition of two CT images using two different acceleration voltages) that provides more constraints on the differentiation of heavy elements in samples. The validity of the proposed method was evaluated theoretically by analyzing an exact solution of the cupping artifact for a homogeneous cylindrical sample of an aqueous solution/suspension containing a heavy element.

Principle of Heavy Element Quantification Using Beam Hardening
In this section, we describe the principle for determining the concentration and atomic number of a heavy element in a sample through quantitative analysis of an X-ray CT image containing a beam hardening artifact. Consider a cylindrical homogeneous sample of 43.0 mM of an iodine aqueous solution (4-cm radius). The linear attenuation coefficient (LAC) spectrum μ of the sample shows a strong dependence on the X-ray photon energy E. This phenomenon is largely caused by the heavy element iodine, which is a strong X-ray absorber (Fig. 1). When this sample is imaged by an X-ray CT system using X-rays accelerated at 100 kV (Fig. 1), a beam hardening artifact (i.e., cupping artifact) appears because of the polychromaticity of the X-ray source spectrum. The cupping artifact in the reconstructed CT image of the homogeneous sample is characterized by a large LAC value near the sample rim and small LAC value at the sample center (Fig. 2a) [13][14][15][16][17]. The reconstructed LAC value distribution in the radial direction (Fig. 2b) is a consequence of the physical interaction of (i) the polychromatic source spectrum i 0 (E) of the X-rays irradiated from the tube and (ii) LAC spectrum μ(E) derived from the sample chemistry. Therefore, if the X-ray source spectrum i 0 (E) is known, it is possible to estimate the sample chemistry (the atomic number and molar concentration of the heavy element) by carefully analyzing the cupping profile in Fig. 2b. This is the main principle of heavy element quantification based on CT image analysis [9,10]. This principle has been applied to a CT image of a 550-mM CeCl 3 aqueous solution in a container [9]. The CT image was obtained experimentally by using a single acceleration voltage of 100 kV. Through numerical simulations of CT image reconstruction obtained by considering various heavy elements for the container image, a trial-and-error search was performed for the atomic numbers of the heavy element candidates that most reasonably reproduced the experimentally obtained cupping profile. This search for appropriate atomic numbers was largely successful. However, analysis using the single-exposure CT image method failed to eliminate heavy element candidates with abnormally high atomic numbers (i.e., 75 Re and 76 Os) relative to the ground truth (i.e., 58 Ce) [9]. This is because the cupping profiles for 75 Re and 76 Os are nearly the same as that for 58 Ce when a CT image is acquired at 100 kV.
In this paper, we propose an improved method to eliminate such candidates through the use of double-exposure CT. This method requires the acquisition of two CT images using two different acceleration voltages. The cupping profile (e.g., Fig. 2b) is dependent on the X-ray source spectrum, as well as the sample chemistry. Therefore, CT images obtained with different acceleration voltages for the X-ray tube provide different constraints that are useful for the identification of an unknown heavy element in a sample. The verification of this concept based on theoretical calculations is presented in the following section.

Heavy Element Quantification Using Double-Exposure CT
In this section, we discuss the determination of the atomic number and concentration of a heavy element in a homogeneous cylindrical sample using double-exposure CT. Determination is performed by systematically calculating the effects of sample chemistry and acceleration voltage on the cupping profile. Cupping profiles that are numerically obtained using image-reconstruction computer software [9,13,17] inevitably contain errors that are introduced by the discretization of sample shapes [13,18] (see the noisy fluctuation of the profile in Fig. 2b). Therefore, an analytical solution [19] was employed for cupping profile analysis in this study to avoid such undesirable errors.  Fig. 1 Representative examples of the photon-energy spectra for an X-ray source accelerated at 100 and 120 kV [14]. Each sum is normalized for unification. The LAC spectra μ(E) for I-H 2 O, Xe-H 2 O, and Ta-H 2 O systems with specific molar concentrations are superimposed. The K absorption edges are indicated by arrows (5.3, 5.5, and 10.8 fJ for 53 I, 54 Xe, and 73 Ta, respectively)

3
The theoretical cupping profile f(r) for a homogeneous cylindrical sample is given as where R is the sample radius, r is the coordinate in the radial direction (Fig. 2a), N is an integer that represents the upper limit of the summation, and Γ is the gamma function [19]. No bowtie filter [20] for beam hardening reduction is assumed. The coefficient C n can be calculated explicitly based on the X-ray source spectrum and LAC spectrum of the sample. For example, the coefficient C 1 is where i 0 (E) is the normalized X-ray source spectrum [19], and E max is the maximum value of the X-ray photon energy. Theoretically E max is determined by the acceleration voltage of the X-ray tube. For example, E max is 100 keV = 16.0 fJ for 100 kV and is 120 keV = 19.2 fJ for 120 kV. The coefficient C n for n ≥ 2 can be readily calculated using a simple recurrence relation with the nth moment of the LAC spectrum of the sample [19].
The continuous distribution f(r) can be reasonably characterized using only two values, namely the reconstructed voxel values at the center [i.e. f(0)] and in the vicinity of the rim of the cylindrical sample [9]. Based on Eq. (1), the limiting behavior at the sample rim [19] is that f(r) → C 1 as r → R. Therefore, it is reasonable to assume that f(r) ≈ C 1 in the immediate vicinity of the sample rim (Fig. 2b).
The values of f(0) and C 1 are thus utilized as measures to characterize the severity of cupping artifacts.
We consider the following simple object as a sample. A simple substance containing a heavy element and pure water is confined in a cylindrical plastic container. The X-ray absorption of the very thin plastic container is negligible compared to that of the heavy element. The sample is a homogeneous aqueous solution or fine-grained suspension, meaning the ion/grain size of the simple substance is very small compared to the voxel dimension and it can be accurately approximated as a homogeneous substance. The effects of sample chemistry and acceleration voltage on the value of C 1 value were calculated systematically for a fixed value of f(0) according to the procedure presented by Nakano and Nakashima [19]. The detailed input data for the exact solutions were as follows. Heavy elements ranging from 47 Ag to 84 Po were examined as simple substances. The radius of the sample R was 4 cm, and N was 100. X-ray source spectra were derived for acceleration voltages of 100 and 120 kV (Fig. 1). The value of f(0) was tentatively fixed at 0.2800 (≈ 200 Hounsfield units (HU)) and 0.2571 cm −1 for the acceleration voltages of 100 and 120 kV, respectively. These parameters correspond to an I-H 2 O mixture (43.0 mM iodine). We followed the method presented by Nakashima and Nakano [9,10] to estimate the bulk density of the mixture, which was required for LAC spectrum calculations based on the National Institute of Standards and Technology database [21]. Three examples of calculated LAC spectra are presented in Fig. 1.
For the iodine-bearing samples, in addition to the line profile for the I-H 2 O mixture, a line profile for a KI-H 2 O mixture (42.2 mM KI) was also calculated. The reconstructed voxel values of the CT images are very sensitive to the atomic number of the heaviest element in the sample and less sensitive to the atomic numbers of lighter elements [12,22]. To demonstrate this lesser sensitivity, an iodine compound containing a lighter element (i.e., potassium) was considered. A 42.2-mM KI solution was considered because the corresponding concentration value yields the same f(0) value (i.e., 0.2800 cm −1 ) for an acceleration voltage of 100 kV.

Results
Figure 2b presents a representative example of the theoretical cupping profile of the I-H 2 O mixture in Fig. 1 (43.0-mM iodine mixture). The theoretical curve f(r) is relatively smooth compared to the noisy line profile obtained numerically by using an image-reconstruction computer program. Therefore, the use of the exact solution from Eq. (1) ensures that the severity of cupping artifacts can be assessed with high accuracy.
As for the 38 heavy elements ranging from 47 Ag to 84 Po, the values of C 1 and the molar concentrations determined using Eq. (1) are plotted in Fig. 3 against the atomic numbers of each heavy element in the aqueous mixture. The data points in Fig. 3 were determined by a trial-anderror search for the C 1 value and molar concentration of each heavy element that most reasonably reproduced the f(0) value of 0.2800 cm −1 for 100 kV. The trial-anderror search was also performed for the f(0) value of 0.2571 cm −1 with respect to 120 kV; the results are plotted in Fig. 4, showing C 1 values and molar concentrations significantly different from those in Fig. 3. The results for the three samples whose LAC spectra were presented in Fig. 1 are also included in Figs. 3 and 4. The characteristic C 1 curves with minimum values in Figs. 3 and 4 are essentially identical to the curve presented in the study by Nakashima and Nakano [9] (see Fig. 6 in the literature).
The curves with minimum values are a consequence of the sensitivity of the severity of cupping artifacts to the relative position between the K absorption edge of a heavy element and the spectral peak of an X-ray source [9,16,17]. All data needed for the material differentiation using double-exposure CT are given in Figs. 3 and 4. The data should be arranged for convenience in terms of the determination of the atomic number and molar concentration. Thus, the data points in Figs container are important for the process of heavy element quantification, as discussed in the following section. In contrast to the exact solutions derived from Eq. (1) in this study, experimentally obtained CT images acquired using a CT apparatus contain undesirable artifacts [23] and noise [24]. Measurement errors caused by the noise and artifacts may reduce the accuracy of the determination of atomic numbers and concentrations of heavy elements [9]. Based on the literature [24] on the experimentally obtained standard deviation of the voxel-value fluctuation for real CT imagery, an error level of 7 HU corresponding to 0.0016 cm −1 was tentatively chosen in the present study. The error level of 0.0016 cm −1 was expressed as the radii of circles in Fig. 5. The dataset in Fig. 3 implies that fluctuation in f values (± 0.0016 cm −1 ) results in fluctuation of molar concentrations of approximately ± 0.8 mM, which is indicated by solid lines in Fig. 6.

Discussion
The quantification of a heavy element in a sample based on the analysis of double-exposure CT images can be performed as follows. Iodine that is one of hazardous heavy elements was chosen as the target element in the present study. We assume that the ground truth is the 43.0-mM iodinewater mixture, which yields values of C 1 = 0.3084 cm −1 and f(0) = 0.2800 cm −1 at 100 kV, and C 1 = 0.2818 cm −1 and f(0) = 0.2571 cm −1 at 120 kV. Therefore, the problem can be defined as using Figs. 5 and 6 to search for the atomic number and molar concentration of the heavy element in the cylindrical aqueous sample with R = 4 cm that reproduce the values of C 1 and f(0) mentioned above. Although the LAC spectra are very different between the 45.7-mM tantalum and 43.0 mM-iodine aqueous samples (see Fig. 1), the sample containing 73 Ta yields nearly the same C 1 value (i.e., 0.3081 cm −1 ) as the sample containing 53 I (i.e., 0.3084 cm −1 ) at 100 kV (Fig. 3). This makes it difficult to eliminate 73 Ta from the heavy element candidates for the sample. It is a poor result to fail to distinguish iodine from tantalum because their atomic numbers are very different. A similar difficulty arises at 120 kV (Fig. 4) for the C 1 value of 80 Hg (i.e., 0.2819 cm −1 ), which is nearly the same as that of 53 I (i.e., 0.2818 cm −1 ). This difficulty was encountered in a previous study [9] on heavy element quantification using single-exposure CT.
This difficulty can be avoided by using double-exposure CT based on Fig. 5. An error circle centered on the iodinewater mixture is drawn in Fig. 5 to represent the ground truth. The circle overlaps only with the adjacent circle centered on the xenon-water mixture. No overlap occurs with the circles for tellurium, tantalum, and mercury. Therefore, the heavy elements that reproduce the values of C 1 and f(0) are iodine and xenon. It may be acceptable to fail to eliminate 54 Xe as a candidate because its LAC spectrum is nearly the same as the 53 I spectrum (see Fig. 1). It should be noted that double-exposure CT enables the successful elimination of tantalum and mercury, while single-exposure CT fails to eliminate these candidates. This is a physical consequence of C 1 being sensitive to the acceleration voltage of the X-ray tube for a fixed value of f(0), which allows it to provide the different constraints required for element identification based on CT at different acceleration voltages.
Once the atomic numbers of the candidates are determined, the molar concentrations can be readily estimated based on Fig. 6. According to Fig. 6, the iodine concentration is 43.0 mM at both 100 and 120 kV, whereas the xenon concentration is 41.0 mM at 100 kV and 40.8 mM at 120 kV. The two data points for iodine and xenon in Fig. 6 fall within the error lines of ± 0.8 mM, which confirms that the choice of these two elements is self-consistent and reasonable. It also should be noted that tantalum and mercury are located out of the error lines of ± 0.8 mM, demonstrating that the two candidates are eliminated again in Fig. 6.
The data points for the KI-H 2 O mixture (42.2 mM potassium iodide) agree well with those for the I-H 2 O mixture (43.0 mM iodine). This is because C 1 values are less sensitive to the atomic numbers of lighter elements, such as 19 K, while they are more sensitive to those of the heaviest elements [22]. This ensures that the accuracy of the doubleexposure CT method is not significantly affected by species of light elements in compounds with heavy elements. Figures 5 and 6 demonstrate the advantage of the present study over the conventional double-exposure CT method. In the conventional method, f(r) is fixed to be almost constant within the homogeneous cylindrical sample (i.e., f(0) ≈ C 1 ) because the cupping artifact is suppressed by specific methods [8,[13][14][15]20]. As a result, Fig. 5 is no longer available, and only the constraint on self-consistency (Fig. 6) can be used for the heavy-element determination. According to Fig. 6, as many as ten candidates ( 47 Ag to 56 Ba) fall within the error lines of ± 0.8 mM. In contrast, the combination of Figs. 5 and 6 allows only 53 I and 54 Xe as candidates, demonstrating the advantage of the method of the present study.
The method for CT-based heavy element quantification that was proposed by Nakashima and Nakano [9,10] is unique in that it utilizes the dependence of the atomic number and concentration on the cupping profile. The accuracy of this method can be improved by employing doubleexposure CT, as shown in Figs. 5 and 6. This is analogous to the advantages of two-dimensional nuclear magnetic resonance (NMR) [25] compared to one-dimensional NMR, where using two dimensions significantly improves the ability to differentiate materials. However, a careful comparison between the image analysis methods proposed in this paper and the conventional double-exposure X-ray imagery [7,[26][27][28][29] should be performed in future. Additional directions for future research include the following four subjects.
(i) We have demonstrated successful heavy element identification using the analytical solution in Eq. (1) for a cylindrical phantom sample. One of the next possible research steps is to analyze real samples (e.g., natural ores or polluted soils with heavy elements). However, the shape of the real samples is not always cylindrical. The applicability of the proposed method to a sample with an irregular shape (noncylindrical shape) should be evaluated. Numerical (not analytical) image reconstruction simulations [9,10] for reproducing experimentally obtained beam hardening could be a possible solution. (ii) Real samples often contain two or more heavy elements as compounds (e.g., Hg 2 I 2 ). While Figs. 3, 4, 5, 6 show that effects of a light element in a compound (i.e., potassium in potassium iodide) are negligible, those of heavy elements (e.g., mercury in mercury iodide) would be significant due to the strong X-ray attenuation. Thus, re-calculation of the all curves in Figs. 3, 4, 5, 6 is needed when applying the method of the present study to compounds with two or more heavy elements. (iii) Regarding lanthanoids, the first derivative of the C 1 curve is nearly zero and the C 1 value calculated using Eq. (2) is less sensitive to the atomic number (see the minimum values in Figs. 3 and 4). Therefore, if the ground truth is a lanthanoid, then estimation accuracy would decrease because many candidates are included in the error circles of 0.0016 cm −1 in Fig. 5. This problem could be solved by perform-ing CT at acceleration voltages other than 100 and 120 kV (e.g., 80 and 140 kV). This is because the position of the C 1 minimum value shifts with the acceleration voltage, as confirmed by a comparison between Figs. 3 and 4. Additionally, the reduction of the radius of the error circle by increasing the tube current could be another solution. (iv) This study is based on the physics of beam hardening that is derived from the non-linear relationship between the X-ray projection (logarithmic intensity of the attenuated X-ray) and the sample penetration depth [19]. Such non-linearity also occurs in X-ray radiography (not tomography) of a sample with an irregular shape [28]. Thus, the method of the present study could be potentially extended to the doubleexposure radiography [29], which is an economically more reasonable apparatus than tomography.

Conclusions
The principle for nondestructive heavy element identification by means of double-exposure X-ray CT was presented for analytical chemistry. This method involves quantitative CT image analysis of the beam hardening resulting from the use of a polychromatic X-ray source. By using analytical expressions for the cupping profile of the beam hardening for a homogeneous cylindrical sample, as well as the dependence of the profile on the atomic number and concentration of a heavy element, and on the acceleration voltage of the X-ray tube, heavy elements can be systematically identified. Careful analysis of these dependencies enabled the atomic numbers and molar concentrations of heavy elements to be determined simultaneously within a certain estimation error range.