Mineral quantitative characterization method based on basis material decomposition model by dual-energy computed tomography

Abstract

BACKGROUND:

Dual-energy computed tomography (DECT) can reconstruct electron density ρ_e and effective atomic number Z_eff distribution for material discrimination. Image-domain basis material decomposition (IBMD) method is a widely used DECT method. However, IBMD method cannot be used for mineral identification directly due to limitations of complex basis material determination, beam hardening artifacts, and inherent errors caused by approximate empirical formulas.

OBJECTIVE:

This study proposes an improved IBMD (IIBMD) method to overcome the above limitations.

METHODS:

In IIBMD method, the composition of basis material is optimized to obtain accurate decomposition coefficients, which enables accurate ρ_e and Z_eff distribution. Moreover, the thickness of basis material is optimized to reduce the effect of beam hardening. Furthermore, two formulas in place of empirical formulas are proposed to calculate ρ_e and Z_eff. Finally, a threshold technique is applied to separate different mineral phases.

RESULTS:

Numerical simulations and practical experiments using a photon-counting detector CT system are implemented to verify IIBMD method. Results show that the relative errors of ρ_e and Z_eff for seven common minerals are down to 5%, lower than most of the existing DECT methods for rocks. Reasonable volume fraction results of mineral phases are thus obtained through threshold segmentation.

CONCLUSIONS:

This study demonstrates that the proposed IIBMD method has high practical value in mineralogical identification.

Keywords

Dual-energy computed tomography basis material decomposition mineral quantitative characterization basis material determination

1 Introduction

X-ray computed tomography (XCT) is a non-destructive testing technique that quantifies the internal components of objects based on differences in linear attenuation coefficients among materials [1]. Recently, XCT is widely applied in geology to obtain quantitative information about component, structure, and pore distribution of rock [2 –5]. However, XCT cannot identify minerals with similar linear attenuation coefficients in heterogeneous rocks. Dual-energy computed tomography (DECT) can overcome this limitation by decomposing linear attenuation coefficient measured at two different voltages into electron density ρ_e and effective atomic number Z_eff [6]. Conventional DECT methods can be classified into three categories: image domain method [7], projection domain method [8, 9], and iterative method [10]. Image domain method reconstructs dual-energy images by conventional CT reconstruction methods and extracts quantitative decomposition results from dual-energy images. Image domain method is applicable for mineral characterization in geology due to its efficiency and applicability. Several studies have attempted to apply the DECT image domain method in geology.

The image-domain DECT method in geology was used by Van Geet [11] firstly, who adapted the mathematical formulas previously elaborated in literature [12] to quantify ρ_e and Z_eff of reservoir rocks. Later, two-parameter model [13] was used to decompose ρ_e and Z_eff from dual-energy images of common minerals in geology [14, 15]. Simple mineralogical segmentation was achieved. However, the performance of the above methods depends highly on calibration with many standard materials, which limits the segmentation accuracy. Recently, the stoichiometric calibration method elaborated for medical purposes was applied to identify different mineral phases [16]. Common minerals in geology were distinguished successfully. However, the stoichiometric calibration method failed in the identification of dense samples. Beam hardening artifacts in reconstructed images cannot be removed completely, either.

Image-domain basis material decomposition (IBMD) method is another typical image domain material decomposition method that is widely used to extract ρ_e and Z_eff from DECT images. The precision of IBMD method [8] can be improved by the optimization of basis material without the need for extensive calibrations. However, the application of IBMD method to geological samples faces some challenges. Firstly, the optimization of basis material, which decides the accuracy of IBMD method [8, 17], has rarely been discussed in geology. Secondly, the beam hardening artifact cannot be eliminated completely, which is not desirable in quantitative imaging. Thirdly, there exist inherent errors caused by the use of approximate empirical formulas for ρ_e and Z_eff, especially for materials with high Z_eff [18], such as geological samples.

In this paper, an improved IBMD (IIBMD) method is proposed to obtain accurate quantitative characterization results of geological samples. In IIBMD method, the composition of basis material is determined pixel by pixel to obtain accurate decomposition results. Moreover, beam hardening artifacts are reduced by optimizing the thickness of basis material. In addition, two formulas in place of the empirical formulas are proposed to calculate ρ_e and Z_eff to minimize the inherent error. These two formulas are established based on virtual monochromatic images at optimum energies (1MeV and 10KeV in this paper). After that, a simple thresholding technique is applied to the images of ρ_e and Z_eff to separate different mineral phases. Simulation and practical experiments were performed on a photon-counting detector (PCD)-based DECT system to validate the IIBMD method.

2 Method

2.1 Image domain basis material decomposition method

Image-domain basis material decomposition (IBMD) method is based on the assumption that the linear attenuation coefficient μ of any material can be expressed as a linear combination of μ of two basis materials [8] at given X-ray energy E. Therefore, IBMD model is expressed as: $μ^{E} = μ_{1}^{E} a_{1} + μ_{2}^{E} a_{2}$ (1) where μ ^E represents μ of any material at energy E, represents μ of basis material i (i = 1, 2) at energy E, and a_i represents the volume fractions of basis material i. In the dual-energy situation, μ of any material at two energies is then expressed as: $μ^{L} = μ_{1}^{L} a_{1} + μ_{2}^{L} a_{2}$ (2) $μ^{H} = μ_{1}^{H} a_{1} + μ_{2}^{H} a_{2}$ (3) where L and H represent low energy and high energy respectively, and subscript 1 and 2 represent the two basis materials. Then a₁ and a₂ are obtained by solving formula (2) and formula (3). Electron density ρ_e and effective atomic number Z_eff are obtained by the following two formulas [8]: $ρ_{e} = a_{1} ρ_{e 1} + a_{2} ρ_{e 2}$ (4) $Z_{eff} = {(\frac{a_{1} ρ_{e 1} Z_{1}^{n} + a_{2} ρ_{e 2} Z_{2}^{n}}{a_{1} ρ_{e 1} + a_{2} ρ_{e 2}})}^{1 / n}$ (5) where ρ_e1, ρ_e2, Z₁, and Z₂ are electron densities and atomic numbers of two basis materials, and n is a constant (n lies between 3 and 4). ρ_e and Z_eff for compounds are defined as follows [8 , 20]: $ρ_{e} = 2 ρ [\frac{\sum_{i = 1}^{t} c_{i} Z_{i}}{\sum_{i = 1}^{t} c_{i} A_{i}}]$ (6) $Z_{e f f} = {(\sum_{i = 1}^{t} β_{i} Z_{i}^{m - 1})}^{\frac{1}{m - 1}}$ (7) where ρ is the density of compound in g/cm³, c_i is the number of atoms for the ith element in a molecule of the compound, Z_i is the atomic number of the ith element, A_i is the atomic weight of the ith element, t is the number of element types in the compound, β_i is the fractional number of the electrons of the ith element, and m is a constant between 3 and 5. It is preferred to set m to 3.6 for materials with Z < 6, and 4.1 for materials with Z > 6. When the chemical formula and density for a compound are known, the theoretical values of ρ_e and Z_eff can be obtained by formula (6) and formula (7).

IBMD method works well in medical applications. However, IBMD method might be problematic for geological samples which are denser and more complex than human tissues. Firstly, the determination of basis material for geological samples is a puzzle, which has rarely been investigated. In addition, IBMD method suffers a lot from beam hardening effect, which hampers the quantitative characterization of material [21]. Furthermore, inherent errors caused by empirical formulas are significant, particularly for geological samples with high Z_eff. To overcome the above limitations, IIBMD method is proposed. In IIBMD method, the composition of basis material is determined pixel by pixel to accommodate the non-homogeneous property of rock. Additionally, the thickness of basis material is optimized to remove beam hardening artifacts, which is also carried out pixel by pixel. In addition, two formulas in place of the empirical formulas are proposed to calculate ρ_e and Z_eff. These two formulas are established from virtual monochromatic images at optimum energies to minimize the inherent error for geological samples with high Z_eff.

2.2 Basis material determination

2.2.1 Basis material composition determination

For IBMD method, the determination of basis material is the primary consideration. Basis material is determined based on component analysis [17, 22] in medical application. The determination of basis material is easy to complete since the composition of medical material is fairly known (e.g., brain tissue, water, fat, bone). However, component analysis is no longer applicable to geological samples with complex composition and random distribution of mineral phase. Therefore, the attenuation property of material at the pixel scale is utilized to determine the composition of basis material. The details are as follows.

To quantify the attenuation property of the material in the dual-energy situation, a metric which is called R - value is defined: $R_{j} = \frac{μ_{j}^{H}}{μ_{j}^{L}} = \frac{F (H, Z_{e f f j})}{F (L, Z_{e f f j})} = G (H, L, Z_{e f f j})$ (8) where j represents the jth material. For monochromatic X-ray source, R - value is a constant that is only related to Z_eff once the energies L and H are determined. When the X-ray source is polychromatic, R - value varies with the thickness of material due to the thickness effect [23]. Even so, a great dependence between R - value and Z_eff still exists [23]. However, an overlap of the R - value exists among materials with extremely different Z_eff, which puzzled the determination of basis material. In Fig. 1 (a), for example, even Z_eff of calcium is significantly different from that of titanium and vanadium, R - value of calcium is still similar to R - value of titanium and vanadium.

Fig. 1

(a) R - values of five pure elements vary with the thickness.

To solve the problem above, $μ_{i}^{L}$ is determined as the second metric of substance attenuation characteristics. As shown in formula (9), a model function is proposed to offer an expression for $μ_{i}^{L}$ as a function of the thickness. Then curve R - μ^L of basis material including thickness information is obtained by formula (8) and formula (9). $μ = m_{1} + \frac{m_{2}}{t} + \frac{m_{3}}{t^{2}}$ (9) where μ represents the polychromatic attenuation value, m_i (i = 1, 2, 3) is a constant, and t represents the thickness of the basis material. m_i (i = 1, 2, 3) in formula (9) is obtained by the least-squares fitting of dataset μ - t, as detailed in Section 3.2.

The curves of R - μ^L show excellent separation capability for materials with different attenuation characteristics. As shown in Fig. 1 (b), calcium, titanium and vanadium are well separated. Once curves R - μ^L of different basis materials are obtained, the basis material determination is performed on the plane R - μ^L. To simplify the operation, an enclosed space, which is called the two-dimensional basis material combination property space (2D BMCP space), is formed by connecting curves R - μ^L of any two basis materials (as shown in Fig. 1 (b)).

The proposed basis material determination method consists of three steps: (i) A basis material library devoted to geology is prepared. (ii) A 2D BMCP space library is created from DECT images. Every element in the library is the 2D BMCP space of two materials picked from the basis material library. (iii) Two optimal basis materials are determined for each pixel of DECT images by searching for one optimal 2D BMCP space in 2D BMCP space library. The optimal 2D BMCP space is the minimum one that can enclose point (μ^L, R - value) of the pixel. If no 2D BMCP space in the library can enclose point (μ^L, R - value), the optimal 2D BMCP space is the closest one to the point (μ^L, R - value). As shown in Fig. 1 (b), the points (μ^L, R - value) of hematite and siderite locate in 2D BMCP space 1 and 2 respectively. Therefore, calcium and titanium are selected as the basis materials for siderite, while titanium and vanadium are the ones for hematite.

2.2.2 Basis material thickness optimization

After the compositions of two basis materials are determined for each pixel, IBMD method is performed to decompose two basis materials. The material decomposition may be imperfect due to beam hardening [21]. Cupping artifacts caused by beam hardening are shown in the image of $μ_{i}^{E}$ and transferred to the image of a_i, manifesting as attenuation values varying with position for a single material rather than a single value for a single material [24]. Many methods to attenuate beam hardening have been proposed, ranging from hardware methods [25] to software methods [8 , 26–28], from the projection domain method [8, 26] to the image domain method [27]. However, IIBMD method is devoted to attenuate beam hardening of rock samples in the process of material decomposition in image domain, which does not require spectral information [8], the prior information of material [28], or curve fitting model with extensive amounts of calibration [26].

Taking single component material for example, the above method is illustrated as follows. Recalling formula (1), IBMD model contains three parameters: μ^E, $μ_{i}^{E}$ , and a_i. Cupping artifact exists in the image of μ^E. In conventional IBMD method, the value of $μ_{i}^{E}$ is the same for each pixel of single component sample. Moreover, basis material decomposition is a linear operation. Thus, cupping artifact in the image of μ^E is transferred to the image of a_i. The core assignment in this paper is to correct the image of $μ_{i}^{E}$ to fit the extent of beam hardening of sample. In this way, “beam hardening artifact” appears in both the image of μ^E and the image of $μ_{i}^{E}$ . The artifact in the image of a_i caused by beam hardening is then corrected after the linear operation.

Corrected image of $μ_{i}^{E}$ is obtained by optimizing the thickness of basis material. Different from the conventional method, $μ_{i}^{E}$ in IIBMD method is the reconstructed value in the center of the phantom which varies with the thickness [23] instead of the average reconstructed value. However, IBMD model alone cannot provide sufficient constraints for determining the thickness of basis material. Therefore, volume conservation constraint is added as another condition. According to volume conservation constraint, two volume fractions satisfy the following two equations: $a_{1} + a_{2} = 1$ (10) $0 < a_{1} < 1, 0 < a_{2} < 1$ (11)

The principle for determining the thickness of basis material is then expressed as follows: the values of $μ_{i}^{E}$ and μ^E satisfy the IBMD model and volume conservation constraint best for each pixel of the images. Figure 2 (a) shows the energy map of siderite sample and basis material that satisfies both IBMD model and volume conservation constraint. In the energy map, cupping artifact appears as scattered point clouds. As shown in Fig. 2 (a), even if the cupping artifact exists, similar a_i is obtained by choosing different thicknesses of basis material (or we can say different $μ_{i}^{E}$ ). Combined with the 2D BMCP space, the thickness of basis material is determined as follows:

Fig. 2

(a) The energy map of siderite and the linear combination of basis material titanium and vanadium with varying thickness t = 2 ∼30mm.

By setting the thickness of basis material as z-axis, 2D BMCP space is extended to three-dimensional (3D) BMCP space (see the 3D surface in Fig. 2 (b)).

The thickness of basis material is determined by mapping 2D points (μ^L, R - value) of the scanned sample onto 3D BMCP space. The optimum thickness is the z-coordinate of each mapping point. Figure 2 (b) provides a graphical illustration for basis material thickness determination. The points (μ^L, R - value) of hematite are mapped to 3D BMCP space of titanium and vanadium to determine the optimum thickness.

2.3 Electron density and effective atomic number calculation

After two basis materials are determined, the images of a₁ and a₂ without the cupping artifact are obtained. Then the next step is to calculate ρ_e and Z_eff. As mentioned above, formula (4) and formula (5) are two empirical formulas that lead to inherent errors in the images of ρ_e and Z_eff. The detailed analysis is shown in literature [18]. To minimize the inherent errors, two formulas based on virtual monochromatic image are established to calculate ρ_e and Z_eff.

Once a₁ and a₂ are determined, the virtual monochromatic image at any energy is given by formula (1). X-ray attenuation for E < 1MeV arises from three physical mechanisms: Coherent scattering, Compton scattering, and Photoelectric absorption [29]. Therefore, μ^E is expressed as follows: $μ^{E} = ρ_{e} (Z_{eff}^{m} f_{Coh} (E) + f_{Compton} (E) + Z_{eff}^{n - 1} f_{Photo} (E))$ (12)

In formula (12), f_Coh (E), f_Compton (E), and f_Photo (E) denote the Coherent scattering coefficient, the Compton scattering coefficient, and the Photoelectric coefficient, which are only related to the energy E, whereas m and n are constants.

Figure 3 shows the variation of μ^E for three different minerals at the range of 0 . 01MeV < E < 1MeV which is calculated according to the NIST XCOM database [30]. At E = 1MeV, as shown in Fig. 3, Compton scattering interaction accounts for more than 99.9% of photon-matter interaction. Therefore, the virtual monochromatic image at E = 1MeV provides a linear relationship between μ^E and ρ_e of the minerals:

$μ^{E_{1}} = ρ_{e} f_{Compton} (E_{1})$ (13) where E₁ = 1MeV. A dataset μ^E₁ - ρ_e is fitted to formula (13), then the energy-dependent coefficient f_Compton (E₁) is obtained.

Thus, ρ_e is derived from formula (13). Similarly, at E = 10KeV, the photon-matter interaction exclusively involves Compton scattering interaction and Photoelectric absorption interaction. Then the attenuation coefficient is approximately expressed as: $μ^{E_{2}} = ρ_{e} (f_{Compton} (E_{2}) + Z_{eff}^{n - 1} f_{Photo} (E_{2}))$ (14) where E₂ = 10KeV. To describe the power-law dependence between μ^E₂ and Z_eff, an exponent U is defined as follows: $U = \frac{μ^{E_{2}}}{μ^{E_{1}}} = k (E_{2} {, E}_{1}) Z_{eff}^{n - 1} + g (E_{2} {, E}_{1})$ (15) where $k (E_{2}, E_{1}) = \frac{f_{Photo} (E_{2})}{f_{Compton} (E_{1})}$ and $g (E_{2}, E_{1}) = \frac{f_{c o m p t o n} (E_{2})}{f_{c o m p t o n} (E_{1})}$ . In formula (15), k (E₂, E₁) and g (E₂, E₁) are determined by fitting dataset U - Z_eff to formula (15). Then Z_eff can be obtained by formula (15).

Fig. 3

(a) (c) Example of μ^E originating from three photon interactions at the range of 0 . 01MeV < E < 1MeV for three diagenetic minerals. (d) (f) Example of the percentage of three photon interactions at the range of 0 . 01MeV < E < 1MeV for three diagenetic minerals.

2.4 Quantitative material characterization method

The workflow of our proposed method is shown in Fig. 4. First, a basis material library devoted to geological samples is developed (See Section 3.1 for detailed information). Every library member is set as a phantom with adjustable thickness. The 2D and 3D BMCP space library is built using the dual-energy imaging data of these basis material phantoms. Then the composition of basis material is determined using 2D BMCP space. After that, 3D BMCP space library is used to optimize the thickness of basis material. ρ_e and Z_eff of the scanned sample are obtained with the procedure in section 2.3.

Fig. 4

The diagram of the proposed method.

The final goal is to derive quantitative mineralogical results from the images of ρ_e and Z_eff. For this purpose, a threshold technique is carried out on the images of ρ_e and Z_eff. The upper and lower thresholds of ρ_e and Z_eff are determined based on the truth value of ρ_e and Z_eff of the minerals (upper/lower threshold = truth value ± truth value × 5 % in this paper). Possible mineral phases are listed in advance based on the knowledge of the geological samples. Pixels with the value of ρ_e and Z_eff between the same upper and lower thresholds are labeled as the same corresponding mineral. Pixels with the value of ρ_e and Z_eff above chalcopyrite were classified as unknown high-density mineral and those with the value of ρ_e and Z_eff below quartz are classified as porous. Thus, through the combination of the information of ρ_e and Z_eff, the material quantitative characterization of the geological sample is completed.

3 Numerical simulation study

3.1 Experiment settings

Numerical simulations were implemented to investigate the performance of the proposed method. The X-ray spectrum was obtained from Spectrum GUI, an open-source energy spectrum generation program. No filter for the CT scanner was selected. Two narrow spectra: 30KeV ∼ 50KeV for the low-energy spectrum and 70KeV ∼ 90KeV for the high-energy spectrum were used for DECT imaging. In the numerical simulation study, a simulation phantom (Fig. 5) consisting of five minerals was imaged. ρ_e and Z_eff for the five minerals above cover a wide range of ρ_e and Z_eff of common minerals. The image size of the phantom is 512 × 512 pixels and each pixel is 100μ m × 100μ m. A basis material library devoted to geological samples was developed. Seven pure chemical materials: 20% sodium chloride solution (20% NaCl), aluminum (Al), silicon (Si), calcium (Ca), titanium (Ti), vanadium (V), and iron (Fe) were selected as the library members according to our prior knowledge. Every library member was set as a model with adjustable thickness and imaged to obtain the dataset μ - t. The attenuation coefficients of the materials used in this study were generated from the NIST XCOM database [25]. DECT projections under two spectra were simulated using 2D fan-beam CT geometry. FBP algorithm was applied to reconstruct the attenuation value distribution from the projections.

Fig. 5

The simulation phantom.

3.2 Simulation results

The unknown coefficient m_i (i = 1, 2, 3) in formula (9) was determined by fitting the dataset μ - t of basis material to formula (9). As shown in Fig. 6, there is an excellent agreement between formula (9) and the simulation data of iron and silicon with the thickness larger than 1 mm. The determination coefficients R² of the fitting for the above six basis materials are all greater than 0.99, which demonstrates the validity of our proposed model in formula (9).

Fig. 6

Simulated dataset μ - t and fitted lines for formula (9) of iron (Fe) and silicon (Si) at low and high energy.

The unknown coefficients: f_Compton (E₁), k (E₂, E₁), and g (E₂, E₁) in formula (13) and formula (15) were determined by fitting the data of several pure chemical elements to formula (13) and formula (15). The detailed data is listed in Table 1. Figure 7 shows the fitting results. It is noted that the determination coefficients are all above 0.999 (R² > 0 . 999), which demonstrates the validity of the proposed formulas.

Table 1

The reference data of the pure chemical elements used in the least squares fit

Material	μ^E₁	μ^E₂	U	ρ _e	Z _eff
C	0.0124	2.7370	43.0069	1.948	6
Na	0.0059	15.5700	255.0213198	0.926	11
Mg	0.0110	21.0600	334.3467153	1.718	12
AL	0.0166	26.2100	426.5641953	2.602	13
Si	0.0148	33.8800	532.6612686	2.323	14
K	0.0054	6.8158	1271.610821	0.808	19
Ca	0.0099	14.4786	1462.479798	1.547	20
Ti	0.0266	110.7000	1878.772881	4.142	22
Fe	0.0471	170.7000	2847.41511	7.319	26

Fig. 7

(a) Dataset μ^E₁ - ρ_e and fitted lines for formula (13); (b) Dataset U - Z_eff and fitted lines for formula (15).

The determination results of basis material composition and thickness are shown in Fig. 8 (a) and 8 (b). The slice images of ρ_e and Z_eff obtained by IIBMD method are shown in Fig. 8 (c) and 8 (d). The performance of IIBMD method and traditional IBMD method in the calculation of ρ_e and Z_eff is compared in Fig. 8 (e) and 8 (f). In the traditional IBMD method, only one basis material combination (silicon and titanium) with uniform thickness was used for material decomposition, while ρ_e and Z_eff were calculated using formula (4) and formula (5). To fully verify the improvement of IIBMD method, the results of ρ_e and Z_eff from IBMD method with basis material composition determination (IBMDC) and IBMD method with basis material composition and thickness determination (IBMDT) are also shown in Fig. 8 (e) and 8 (f). Obviously, unacceptable errors are shown in the results from IBMD method (yellow lines in Fig. 8 (e) and 8 (f)), which is reduced to an acceptable range by basis material composition determination (blue lines in Fig. 8 (e) and 8 (f)). In addition, the determination of the thickness for basis material significantly attenuates the errors due to beam hardening (blue and green lines in Fig. 8 (e) and 8 (f)). Finally, the results of ρ_e and Z_eff obtained by IIBMD method are closer to the theoretical values than those of IBMDT method (green and red lines in Fig. 8 (e) and 8 (f)). Therefore, formula (13) and formula (15) are effective in improving the calculation accuracy of ρ_e and Z_eff.

Fig. 8

(a) The determination results of basis material composition;

After ρ_e and Z_eff were obtained, the segmented image was obtained by applying threshold technique, as shown in Fig. 9. Table 2 summarizes the results of ρ_e, Z_eff, and volume fraction. The experimental values show an excellent agreement with the theoretical ones within a relative error of 5%, which validates the effectiveness of the proposed method.

Fig. 9

Segmented image of the simulation phantom.

Table 2

Comparison of experimental value and theoretical value for ρ_e, Z_eff, and V of the simulation phantom in Fig. 5

Mineral	ρ _e			Z _eff			V
	Theoretical value	Experimental value	Error (%)	Theoretical value	Experimental value	Error (%)	Theoretical value	Experimental value	Error (%)
Quartz	2.646	2.644	–0.01	11.75	11.87	0.99	87.75	87.46	–0.33
Dolomite	2.824	2.800	–0.85	13.64	13.04	–4.44	3.61	3.68	1.94
Orthoclase	2.549	2.571	0.84	13.34	12.98	–2.64	3.69	3.71	0.54
Corundum	3.433	3.333	–2.89	11.28	11.57	2.65	3.63	3.55	2.20
Hematite	4.988	4.894	–1.89	23.38	23.86	2.05	1.32	1.32	<0.1

*V means the volume fraction and the theoretical values of ρ_e and Z_eff are obtained by formula (6) and formula (7).

4 Practical experiments

4.1 Configuration of the PCD-based DECT system

All the practical experiments were performed on a self-built PCD-based DECT system (Fig. 10), which is developed by The State Key Laboratory of Precision Measuring Technology and Instruments of Tianjin University in China. The DECT system is composed of an X-ray source (Hamamatsu L8121-52), a rotation stage, a CdTe semiconductor PCD (XC-Actaeon FX3), and the controlling system. The detector has 256 × 256 pixels with a pixel size of 100μ m × 100μ m. By setting different energy thresholds, PCD can count the photons above different energy thresholds to obtain CT images of different energy bins. In dual-energy mode, the detector could set two energy thresholds in one acquisition. The low energy threshold amounts to the lowest threshold for the detector to count the photons of the total energy image. The high energy threshold amounts to the lowest threshold for the detector to count the photons of the high energy image. Subtraction operation is applied to obtain the counts of low energy image.

Fig. 10

The PCD-based DECT system.

X-ray source was operated at 55kV and 115mA. Two thresholds used for acquisition were 10kV and 27kV. Thus, two narrow energy bins: 10kV ∼ 27kV for low-energy and 27kV ∼ 55kV for high-energy were achieved. The main dual-energy scanning parameters are summarized in Table 3.

Table 3

Main dual-energy scanning parameters

Parameter	Value
Tube voltage	55kV
Tube current	150mA
Additional filter	None
Source to detector distance	210mm
source to rotation axis distance	125mm
Acquisition frequency	5Hz
Detector resolution	256 × 256 with 100μm × 100μm pixels
Detector material	2.0mm CdTe

4.2 Basis material phantom and geological sample

Six basis materials: graphite, 20% NaCl, aluminum (Al), silicon (Si), titanium (Ti), and iron (Fe) were used for basis material decomposition in the practical experiment. These basis material phantoms are all cylinders of 99.9% purity elementary substances except 20% NaCl. The diameters of these cylinders are 4 mm, 8 mm, 10 mm, 12 mm, and 14 mm. 20% NaCl was placed in centrifuge tubes with diameters of 6 mm, 8 mm, 10 mm, 11 mm, and 13 mm. Seven mineral samples (Table 4, columns 1) and two rock samples (Fig. 11) were used for the validation of IIBMD method. The first rock sample is an artificial sample consisting of quartz, calcite, and siderite. The second one is a natural igneous sample. The possible mineral phases in igneous are classified into six categories (see Table 6 for corresponding parameters) according to prior knowledge [31].

Table 4
Comparison of experimental value and theoretical value for ρ_e and Z_eff of seven mineral samples

Mineral ρ _e Z _eff

Theoretical value Experimental value Error (%) Theoretical value Experimental value Error (%)

Orthoclase 2.548 2.480 –2.84 13.34 12.93 –3.08

Quartz 2.646 2.560 –3.40 11.75 12.19 3.72

Calcite 2.707 2.721 0.48 15.62 15.77 0.96

Dolomite 2.824 3.020 6.94 13.64 13.57 –0.51

Corundum 3.433 3.042 –11.38 11.28 12.02 6.52

Siderite 3.828 3.678 –3.93 20.95 21.55 2.88

Pyrite 4.738 4.678 –1.27 21.95 23.13 5.38

Mineral	ρ _e	Z _eff
Orthoclase	2.548	2.480	–2.84	13.34	12.93	–3.08
Quartz	2.646	2.560	–3.40	11.75	12.19	3.72
Calcite	2.707	2.721	0.48	15.62	15.77	0.96
Dolomite	2.824	3.020	6.94	13.64	13.57	–0.51
Corundum	3.433	3.042	–11.38	11.28	12.02	6.52
Siderite	3.828	3.678	–3.93	20.95	21.55	2.88
Pyrite	4.738	4.678	–1.27	21.95	23.13	5.38

*The theoretical values of ρ_e and Z_eff are obtained by formula (6) and formula (7).

Table 5

Comparison of the theoretical values and experimental values of volume fractions of different phases in the artificial sample

mineral	Theoretical value (%)	Experimental value (%)	Error (%)
Quartz	29	34.26	18.14
Calcite	29	29.71	2.44
Siderite	42	33.69	–19.78
Others	0	2.34	/

Fig. 11

(a) Basis materials; (b) Geological samples: seven mineral samples and two rock samples.

Table 6

Minerals phases in igneous sample with their chemical composition and theoretical ρ_e and Z_eff

Mineral	Chemical composition	ρ _e	Z _eff
Quartz	SiO₂	2.646	11.75
Orthoclase	K[AlSi₃O₈]	2.549	13.34
Andesine	(CaAl₂Si₂O₈)₄₀(NaAlSi₃O₈)₆₀	2.648	13.25
Labradorite	(CaAl₂Si₂O₈)₆₀(NaAlSi₃O₈)₄₀	2.678	16.61
Augite	(Ca,Mg,Fe2+,Fe3+,Ti,Al)₂(Si,Al)₂O₆	3.36	14.85
Magnetite	Fe₃O₄	4.55	21.32

4.3 Practical experiment results

The coefficients m_i (i = 1, 2, 3) in formula (9) were calibrated by fitting dataset μ - t of basis material phantom. The fit for each basis material phantom in low and high energy presents R² > 0.99. ρ_e and Z_eff of seven mineral samples were calculated using IIBMD method. The relative difference between the experimental and theoretical values of ρ_e and Z_eff is shown in Table 4. The experimental values agree well with the theoretical values within one standard deviation, indicating that the IIBMD method is also valid in real situations.

Figure 12 shows the performance of IIBMD method and traditional IBMD method in the calculation of ρ_e and Z_eff for the artificial sample. In the traditional IBMD method, only one basis material combination (silicon and titanium) with uniform thickness was used for material decomposition, while ρ_e and Z_eff were calculated using formula (4) and formula (5). In Fig. 12, the results of ρ_e and Z_eff from IBMD method with basis material composition determination (IBMDC) and IBMD method with basis material composition and thickness determination (IBMDT) are compared to those from IIBMD method. Obviously, unacceptable errors are shown in the results from IBMD method (yellow lines in Fig. 12). The errors are decreased to an acceptable range in IBMDC method (blue lines in Fig. 12), which illustrates the important role of basis material composition determination. In addition, basis material thickness determination significantly attenuates the error due to beam hardening (blue and green lines in Fig. 12). Finally, the results of ρ_e and Z_eff obtained by IIBMD method agree better with the theoretical values than those obtained by IBMDT method. Therefore, formula (13) and formula (15) are effective in improving the calculation accuracy of ρ_e and Z_eff.

Fig. 12

(a) The slice image of ρ_e from IIBMD method; (b) The slice image of Z_eff from IIBMD method; (c) Horizontal profile of ρ_e; (d) Horizontal profile of Z_eff; (e) Vertical profile of ρ_e; (f) Vertical profile of Z_eff.

The spatial distributions (Fig. 13 (b)) and volume fractions (Table 5) of dolomite, calcite, and siderite in the artificial sample were obtained by applying a threshold technique to image ρ_e and image Z_eff. The theoretical volume fractions of different phases were obtained by the drainage method. Although the three mineral phases were separated successfully as shown in Fig. 13 (b), the accuracy of the volume fraction needs to be improved. Figure 14 shows the segmentation result of an igneous sample obtained by the proposed method. The volume fractions of all possible mineral phases are listed in Fig. 14 (b). The experimental value is quite close to the theoretical value, which indicates a reliable segmentation result. Note that the theoretical value in Fig. 14 (b) is obtained by referring to literature [31].

Fig. 13

(a) Low-energy reconstruction image of the artificial sample; (b) 3D segmentation result of the artificial sample.

Fig. 14

3D segmented result of igneous sample. (a) 3D segmented image (see Fig. 14 (b) for corresponding colors). (b) Volume fractions of the separated mineral phases. (c)-(j) Separated image of different phases: (c) quartz; (d) orthoclase; (e) andesine; (f) labradorite; (g) augite; (h) magnetite; (i) other phase; (j) pore and cracks.

5 Discussion and conclusion

The IIBMD method gives accurate quantitative information about ρ_e and Z_eff by automatically determining the optimal basis material and separates different mineral phases using a simple threshold technique. Compared with other methods [14 –16], the IIBMD method does not require the calibration of large number of standard materials. In addition, the IIBMD method also attenuates the effects of beam hardening and performs well on dense samples.

Simulation and practical experimental results verify the potential of the IIBMD method for quantitative characterization of mineral phases. The fitting results of formula (9), formula (13), and formula (15) achieved high relevant coefficients, which indicates that the model functions proposed are accurate. In the simulation experiment, the errors of ρ_e and Z_eff between the experimental and theoretical values are less than 5%. The rock phantom image is accurately segmented into corresponding mineral phases. In the practical experiment, five minerals have a relative error of less than 5% between the experimental and theoretical value of ρ_e and Z_eff. However, a wider difference between the theoretical and experimental value of ρ_e and Z_eff appear in corundum. Several reasons can explain as follows.

First, the main reason is that the optimal basis material is not provided in the practical experiment (aluminum and calcium for Corundum). Second, the energy dependence of the decomposition coefficients also leads to errors. The IIBMD method is based on the premise that a₁ and a₂ in formula (1) are energy-independent. In practice, however, a₁ and a₂ vary with energy. Thus, the synthesized values of attenuation coefficients using the IBMD model are expected to suffer from energy-dependent systematic errors. A detailed analysis is shown in the literature [32]. Even so, the above problems can be solved by continuing to optimize two basis materials. The energy dependence of a₁ and a₂ has been proved to be neglected when suitable basis material is determined.

Moreover, the volume fractions of the artificial sample show a relative error of up to 19% (Table 5). Several reasons could explain the errors. First, gelatinous substance stuck to the surface of the fabricated sample (as shown by the red arrows in Fig. 13) caused misidentification. Therefore, the volume fractions of quartz and calcite are larger than expected. Second, natural minerals contain impurities and pores, which easily affect the segmentation results. For the above problems, optimizing the basis material library to obtain more accurate ρ_e and Z_eff or improving the segmentation algorithm may help.

In summary, IIBMD method demonstrates the potential in 3D quantitative characterization of the structure and components of geological samples. The current study is only a simple validation of IIBMD method. Our future work will devote to optimizing the basis material library and scanning energy to improve the quality of the IIBMD method for common minerals. In addition, an improved segmentation algorithm considering pore and micro porosity (pore whose size is smaller than XCT resolution) will also be included in future quantitative characterization procedures.

Footnotes

Acknowledgments

This work was support by the National Natural Science Foundation of China (NNSFC) [No. 61771328] and the National Key Research and Development Program of China [No. 2017YFB1103900].

References

Cnudde

, Masschaele

, Dierick

, et al., Recent progress in X-ray CT as a geosciences tool, Applied Geochemistry 21 (2006), 826–832.

Siddiqui

, Khamees

A.A.

, Dual-energy CT-scanning applications in rock characterization, in: SPE annual technical conference and exhibition, OnePetro, 2004.

Victor

R.A.

, Prodanović

, Torres-Verdín

, Monte Carlo approach for estimating density and atomic number from dual-energy computed tomography images of carbonate rocks, Journal of Geophysical Research: Solid Earth 122 (2017), 9804–9824.

Cnudde

, Boone

M.N.

, High-resolution X-ray computed tomography in geosciences: A review of the current technology and applications, Earth-Science Reviews 123 (2013), 1–17.

Ketcham

R.A.

, Carlson

W.D.

, Acquisition, optimization and interpretation of X-ray computed tomographic imagery: applications to the geosciences, Earth-Science Reviews 27 (2001), 381–400.

Engler

, Friedman

W.D.

, Review of dual-energy computed tomography techniques, Materials Evaluation 48 (1990), 623–629.

Maaß

, Baer

and Kachelrieß

, Image-based dual energy CT using optimized precorrection functions: A practical new approach of material decomposition in image domain, Medical Physics 36 (2009), 3818–3829.

Zhang

, Cheng

, Zhang

, Chen

, Xing

, A practical reconstruction method for dual energy computed tomography, Journal of X-Ray Science and Technology 16 (2008), 67–88.

Cong

, Man

B.D.

, Wang

, Projection decomposition via univariate optimization for dual-energy CT, Journal of X-Ray Science and Technology 30 (2022), 725–736.

10.

, Zhao

, Zhang

, Accurate iterative FBP reconstruction method for material decomposition of dual energy CT, IEEE Transactions on Medical Imaging 38 (2018), 802–812.

11.

Van Geet

, Swennen

and Wevers

, Quantitative analysis of reservoir rocks by microfocus X-ray computerised tomography, Sedimentary Geology 132 (2000), 25–36.

12.

Coenen

, Maas

, Material classification by dual-energy computerized X-ray tomography, in: Proc. Int. Symp. on Computerized Tomography for Industrial Applications, 1994, pp. 8–10.

13.

Walter

D.J.

, Tkaczyk

E.J.

, Wu

, Accuracy and precision of dual energy CT imaging for the quantification of tissue fat content, in: Medical Imaging 2006: Physics of Medical Imaging, International Society for Optics and Photonics, 2006, pp. 61421G.

14.

Alves

, Lima

, de Assis

, et al., Mineralogy evaluation and segmentation using dual-energy microtomography, X-Ray Spectrometry 44 (2015), 99–104.

15.

Alves

, Lima

, Lopes

, Methodology for attainment of density and effective atomic number through dual energy technique using microtomographic images, Applied Radiation and Isotopes 89 (2014), 6–12.

16.

Martini

, Francus

, Di Schiavi Trotta

and Després

, Identification of Common Minerals Using Stoichiometric Calibration Method for Dual-Energy CT, Geochemistry, Geophysics, Geosystems 22 (2021), e2021GC009885.

17.

Tang

, Ren

, On the conditioning of basis materials and its impact on multimaterial decomposition-based spectral imaging in photon-counting CT, Medical Physics 48 (2021), 1100–1116.

18.

Zheng

, Hao

, Xing

, Methodic error analysis of basis material decomposition method in dual-energy computed tomography, CT Theory and Applications 20 (2011), 153–162.

19.

Chen

, Zhang

, Xia

, et al., Dual-energy CT imaging over non-overlapping, orthogonal arcs of limited-angular ranges, Journal of X-Ray Science and Technology 29 (2021), 975–985.

20.

Vienna, Austria. XMuDat: Photon attenuation data on PC. IAEA-NDS-195. http://www-nds.iaea.org/publications/iaeands/iaea-nds-0195.htm, 1998..

21.

Bam

L.C.

, Miller

J.A.

, Becker

, Basson

I.J.

, X-ray computed tomography: Practical evaluation of beam hardening in iron ore samples, Minerals Engineering 131 (2019), 206–215.

22.

G.H.

, Wang

S.L.

, Zhang

W.B.

, Xu

Z.Y.

, A review of dual energy CT: principles, applications, and future outlook, CT Theory and Applications 21 (2012), 367–386.

23.

Chen

Z.Q.

, Zhao

, Li

, A curve-based material recognition method in MeV dual-energy X-ray imaging system, Nuclear Science Techniques 27 (2016), 1–8.

24.

Brooks

R.A.

, Di Chiro

G.b.

, Beam hardening in X-ray reconstructive tomography, Physics in Medicine 21 (1976), 390.

25.

Meganck

J.A.

, Kozloff

K.M.

, Thornton

M.M.

, et al., Beam hardening artifacts in micro-computed tomography scanning can be reduced by X-ray beam filtration and the resulting images can be used to accurately measure BMD, Bone 45 (2009), 1104–1116.

26.

Herman

G.T.

, Correction for beam hardening in computed tomography, Physics in Medicine & Biology 24 (1979), 81.

27.

Sarkar

, Wahi

and Munshi

, An empirical correction method for beam-hardening artifact in Computerized Tomography (CT) images, J NDTE & International 102 (2019), 104–113

28.

Kyriakou

, Meyer

, Prell

, Kachelrieß

, , Empirical beam hardening correction (EBHC) for CT, J , Medical Physics 37 (2010), 5179–5187.

29.

Jackson

D.F.

, Hawkes

D.J.

, X-ray attenuation coefficients of elements and mixtures, Physics Reports 70 (1981), 169–233.

30.

Hubbell

J.H.

, Berger

M.J.

, Seltzer

S.M.

, et al., Photon Cross Sections Database, https://dx.doi.org/10.18434/T48G6X, November 2010..

31.

Department of Mineralogy and Petrology, Department of Geology, Nanjing University, Igneous Petrology, 1980.07.

32.

Goh

, Liew

, Hasegawa

, Energy-dependent systematic errors in dual-energy X-ray CT, IEEE Transactions on Nuclear Science 44 (1997), 212–217.