Iterative Reconstruction of the Electron Density and Effective Atomic Number using a Non-Linear Forward Model

X-ray computed tomography (CT) is widely used to reconstruct the spatially varying linear attenuation coefficient (LAC) of an object [1]. The spectral response function that is dependent on the X-ray source and the detector is used to quantify the variation in X-ray intensity over the broad range of emitted X-ray energies. Most CT systems have differing spectral response functions that causes the reconstructed LAC to be dependent on the spectral response of the X-ray scanner. Hence, it is not possible to compare the LAC values reconstructed using different X-ray scanners. Also, the estimated LAC does not provide direct information on the density and atomic composition of an object [2].

Dual-energy X-ray computed tomography (DECT) is useful to estimate system independent properties that are independent of the CT system spectral response and other system specifications. Our laboratory based micron scale DECT system scans the same object twice using two different spectra. Using DECT scans, an approach called SIRZ was proposed [3] to estimate the system independent properties of electron density, $\rho_{e}$, and effective atomic number, $Z_{e}$. The electron density, $\rho_{e}$, is a direct measure of the object’s density and the effective atomic number, $Z_{e}$, is useful to predict the atomic composition. In Champley at al.[2], the authors presented the SIRZ-2 approach that improved upon the original SIRZ to produce more accurate estimates for $\rho_{e}$ and $Z_{e}$. SIRZ-2 can be broadly split into four different steps:

1. 1.

   Preprocessing steps that include detector gain correction, detector deblur, and scatter correction.

2. 2.

   Dual-energy basis decomposition to estimate projections of the object at two different synthesized monochromatic basis (SMB) energies.

3. 3.

   Filtered back-projection to reconstruct the LAC values at the two SMB energies,

4. 4.

   Estimation of electron density and effective atomic number from the LACs at the two SMB energies.

The second and fourth steps of SIRZ-2 are based on simplifying approximations that limit its accuracy.

In this paper, we present a new approach called SIRZ-3 to reconstruct the system independent properties of $\rho_{e}$ and $Z_{e}$ while avoiding the approximations made in SIRZ-2. SIRZ-3 replaces the steps 2, 3, and 4 of SIRZ-2 using a single iterative optimization algorithm. Features of the SIRZ-3 method include:

• •

  Non-linear differentiable forward model that expresses the measured DECT data as a direct function of the object’s spatial distribution of $\rho_{e}$ and $Z_{e}$.

• •

  Formulation of an error function that quantifies the deviation of measured data from the output of the forward model given the unknowns, $(\rho_{e},Z_{e})$.

• •

  Method for reconstructing $(\rho_{e},Z_{e})$ from the measured data by minimizing the error function using the limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) solver [4, 5].

The LAC quantifies the magnitude of X-ray attenuation within an object. The LAC varies as a function of the position inside the object, density, atomic composition, and X-ray energy. The dependence of LAC, denoted as $\mu\left(Z,\rho_{e},E\right)$, on the electron density, $\rho_{e}$, and atomic number, $Z$, is given by,

where $E$ is the X-ray energy, $\sigma_{e}(\cdot)$ is the total X-ray electronic cross-section.[2] Equation (1) is only defined for integer valued $Z$ that correspond to pure elements.

In order to represent compounds or mixtures of elements, the cross-section for non-integer $Z_{e}$ is defined as follows [3],

where $Z^{\prime}=\lfloor Z_{e}\rfloor$, and $\delta=\left(Z_{e}-Z^{\prime}\right)$. Note that $Z_{e}$ for a given compound or mixture is not a fundamental physical constant and there are several competing techniques to define $Z_{e}$. We adopt the definition of $Z_{e}$ from Champley et al. [2] since it is a function of the X-ray cross-section with superior system invariance properties. The $Z_{e}$ for a compound or mixture is defined as,

where $\sigma_{c}(E)$ is the tabulated[6] electron cross-section of the material, $S(E)$ is the spectral response at energy $E$, and $M$ is the areal electron density. This definition in (3) is demonstrated to be relatively insensitive to $S(E)$ and the X-ray energy range of interest (30-200keV). For additional details on the definition of $Z_{e}$ and its estimation, the reader may refer to Champley et al. [7].

In this section, we formulate a forward model that mathematically expresses the measured DECT data as a function of the object properties, $\rho_{e}$ and $Z_{e}$. This forward model is necessary to solve the inverse problem of reconstructing the object’s $\rho_{e}$ and $Z_{e}$.

The forward model is derived in discrete coordinate space. Let $\rho_{e,j}$ and $Z_{e,j}$ represent the $\rho_{e}$ and $Z_{e}$ values respectively at the $j^{th}$ voxel inside the object. Then, the LAC at voxel $j$ is given by $\mu_{j,k}=\rho_{e,j}\sigma_{k}\left(Z_{e,j}\right)$, where $\sigma_{k}\left(Z_{e,j}\right)=\sigma_{e}\left(Z_{e,j},E_{k}\right)$ from (2) and $E_{k}$ is the X-ray energy at the $k^{th}$ energy bin. The linear projections (line-integration) of the LAC is expressed as,

where $A_{i,j}$ is the element along the $i^{th}$ row and $j^{th}$ column of the matrix $A$. At the detector, the measurement is the average transmission function $\exp\left(-p_{i,k}\right)$ weighted by the spectral density function $S_{k}$. Hence, the effective transmission is given by,

where $S_{k}$ is the spectral density at the $k^{th}$ energy bin ($S_{k}$ is a discretized version of $S(E)$ used in (3)) such that $\sum_{k}S_{k}=1$. The spectral density function $S_{k}$ quantifies both the spectral density of the X-ray source and the spectral response of the detection system. We compute the negative logarithm of (5), which gives us the measurement data in attenuation space, $-\log\left(t_{i}\right)$.

For DECT, the complete forward model is,

where $S_{k}^{L}$ is the spectral density for X-rays at the low-energy spectrum, $S_{k}^{H}$ is the spectral density at the high-energy spectrum, and $\tilde{y}^{L}_{i}$ and $\tilde{y}^{H}_{i}$ are the forward model predictions for the low-energy and high-energy measurements respectively. To ensure differentiability of the forward model in (6) and (7), we need analytical derivatives for $\sigma_{k}\left(Z_{e,j}\right)$ with respect to $Z_{e,j}$. This derivative is given by,

where $Z^{\prime}_{e,j}=\lfloor Z_{e,j}\rfloor$ and $\left\lfloor\cdot\right\rfloor$ is the floor function.

SIRZ-3 is a framework to solve the problem of reconstructing $\rho_{e,j}$ and $Z_{e,j}$ from measurements $y^{L}_{i}$ and $y^{H}_{i}$ using a forward model that accurately defines an analytical relation between $\left(y^{L}_{i},y^{H}_{i}\right)$ and $\left(\rho_{e,j},Z_{e,j}\right)$ for all $i$ and $j$. SIRZ-3 adopts an iterative approach, where the solution for $\left(\rho_{e,j},Z_{e,j}\right)$ is iteratively refined such that the output of the forward model (whose input is the current estimate for $\rho_{e,j}$ and $Z_{e,j}$) gets progressively closer to the measurements $\left(y^{L}_{i},y^{H}_{i}\right)$. In addition to a forward model, the solution for $\left(\rho_{e,j},Z_{e,j}\right)$ may also be required to satisfy certain sparsity enforcing regularization criteria using a prior model. In this paper, however, we do not use a prior model.

The analytical relation for the forward model used in this paper is shown in equations (6) and (7). The accuracy of SIRZ-3 is dependent on the accuracy of this analytical relation that mathematically describes the measurement process of our X-ray system. Future efforts to improve SIRZ-3 may potentially focus on more accurate forward models that improve upon equations (6) and (7) or more accurate preprocessing algorithms so that (6) and (7) does accurately model the data.

To formulate the reconstruction algorithm, we define the following error function,

where $y^{L}_{i}$ and $y^{H}_{i}$ are the attenuation space measurements at the low-energy and high-energy spectra respectively, and $\tilde{y}_{i}^{L}$ and $\tilde{y}_{i}^{H}$ are the forward model predictions from (6) and (7) respectively. The weights for the penalty terms are given by $w^{L}_{i}=\exp\left(-y^{L}_{i}\right)/N$ and $w^{H}_{i}=\exp\left(-y^{H}_{i}\right)/N$, where $N$ is the number of sinogram pixels summed over all the views. Substituting (6) and (7) in (9), we get,

where $Z_{e}$ is a vector of all $Z_{e,j}$ voxel values, $\rho_{e}$ is a vector of all $\rho_{e,j}$ voxel values. Note that the projection matrix entries $A_{i,j}$ are stored in a memory-efficient sparse representation that only stores the non-zero values and its locations. We used LTT [7] to generate $A_{i,j}$. Then, the reconstruction of SIRZ-3 is obtained by solving,

where $\hat{Z}_{e}$ and $\hat{\rho}_{e}$ are reconstructions of $Z_{e}$ and $\rho_{e}$ respectively, $\rho_{e}^{min}$ and $\rho_{e}^{max}$ are the lower and upper limits for $\rho_{e,j}$ respectively, and $Z_{e}^{min}$ and $Z_{e}^{max}$ are the lower and upper limits for $Z_{e,j}$ respectively. In this paper, we set $Z_{e}^{min}=1$, $Z_{e}^{max}=118$, $\rho_{e}^{min}=0$, and $\rho_{e}^{max}=9.018507\times 10^{-3}$electrons$\times$mol/mm³.

In order to improve conditioning and ensure fast convergence, we optimize z-score normalized versions of $Z_{e}$ and $\rho_{e}$. Let $Z_{e}^{avg}$ and $Z_{e}^{std}$ be the mean and standard deviation of the SIRZ-2 reconstructed voxel values for $Z_{e}$. Similarly, let $\rho_{e}^{avg}$ and $\rho_{e}^{std}$ be the mean and standard deviation of the SIRZ-2 reconstructed voxel values for $\rho_{e}$. Then, we solve,

where we use vector-scalar multiplication¹¹1Vector-scalar multiplication is point wise multiplication of a scalar to every element of a vector. and vector-scalar addition²²2Vector-scalar addition is point wise addition of a scalar to every element of a vector. within the error function $E\left(\cdot\right)$. The SIRZ-3 reconstruction is given by,

We use the python programming language and pytorch framework for implementing and solving (12). Pytorch’s algorithmic differentiation is used to compute the partial derivatives of $E\left(\cdot\right)$ with respect to $Z^{nm}_{e,j}$ and $\rho^{nm}_{e,j}$. We solve equation (12) using the L-BFGS [5, 4] optimization algorithm. Upon convergence and using (13), L-BFGS yields the reconstruction of the electron density, $\hat{\rho}_{e}$, and the effective atomic number, $\hat{Z}_{e}$. We use the LBFGS implementation by Shi et al. [8]. For LBFGS, we use a history size of $64$, weak Wolfe line search, and a convergence criteria that stops the LBFGS iterations when the percentage change in both $Z_{e}$ and $\rho_{e}$ goes below $0.2\%$ for $10$ consecutive iterations. Note that it is important to have a sufficiently strong convergence criteria along with sufficiently large history size to ensure convergence. We use the SIRZ-2 estimates for $Z_{e}$ and $\rho_{e}$ to initialize the optimization in (12).

Using simulated data, we compared the performance of SIRZ-2 and SIRZ-3 for a wide range of materials and object thicknesses. We chose various materials including pure elements, compounds, and mixtures with effective atomic numbers, $Z_{e}$, in the range of $5-30$ since they are among the most commonly occurring materials that are sensitive to X-rays in our chosen energy range. Our choice of object’s thicknesses was such that the maximum attenuation value for the low-energy DECT spectrum was in the range of $0.5-5.5$. Outside this attenuation range, it is typically difficult to extract information from X-ray measurements due to the presence of noise. Such a wide selection of materials and attenuations is useful to demonstrate the trends in performance as a function of $Z_{e}$ and object thicknesses.

We simulated fan-beam X-ray DECT data at $256$ different view angles and $256$ pixel bins. The low-energy X-ray spectrum is generated using a $100$kV bremsstrahlung X-ray source and a $2$mm Aluminum filter. The high-energy X-ray spectrum is generated using a $160$kV bremsstrahlung X-ray source, a $2$mm aluminum filter, and a $2$mm copper filter. DECT data is simulated for a circular object that occupies a $224\times 224$ cross-sectional pixel area within the full field-of-view of $256\times 256$. The source-to-object (SOD) distance and source-to-detector (SDD) distance were chosen to be $500$mm and $1000$mm respectively. The diameter of the object was adjusted such that we obtained maximum attenuation values of $0.5$ to $5.5$ in increments of $1.0$ for the low-energy spectrum. For each attenuation value, given the atomic composition of the object’s material, we computed the diameter that produces the desired low-energy attenuation along the X-ray path passing through the center of the object. Then, the pixel size for the object was obtained by dividing this diameter by $224$ (number of pixels occupied by the object along any direction). We assume a Perkin-Elmer (PE) detector for the detector spectral response. For generating simulated data, we used analytic ray-tracing through geometric solids that is implemented in LTT [7]. We add $0.1\%$ noise to the measurements in transmission space.

For quantitative comparison of $Z_{e}$ and $\rho_{e}$ reconstructions, we utilize two forms of error metrics. The first error metric estimates the bias in the reconstructed values by computing the relative error for either $Z_{e}$ or $\rho_{e}$ as,

where $K\in\left(2,3\right)$ is used to denote either SIRZ-2 or SIRZ-3 methods, $x\in\left(\rho_{e},Z_{e}\right)$ denotes either $\rho_{e}$ or $Z_{e}$, $\mathcal{M}$ is the set of all pixel indices in the interior of the object that excludes the object boundaries, $\left|\mathcal{M}\right|$ is the cardinality (number of elements) of the set $\mathcal{M}$, and $x^{GT}$ is the ground-truth value for $x_{j}$. The second error metric is the relative root mean squared error (RMSE), which is defined as,

The relative RMSE metric combines both the bias and standard deviation of the difference between the reconstruction and the ground-truth.

In Fig. 1, we compare the performance of SIRZ-2 and SIRZ-3 for a disc of copper (Cu) with a diameter of $4mm$ and a maximum attenuation of $4.5$ at the low-energy spectrum. Fig. 1 (a, b) show the ground-truth $\left(Z_{e},\rho_{e}\right)$ images and Fig. 1 (c, d) show the low-energy and high-energy sinograms. SIRZ-2 and SIRZ-3 reconstructions are shown in Fig. 1 (e, f) and Fig. 1 (g, h) respectively. We can see that both the relative error and the relative RMSE reduces with SIRZ-3 when compared to SIRZ-2. By comparing Fig. 1 (e,f) with Fig. 1 (a,b), we see that SIRZ-2 over estimates the $Z_{e}$ and under estimates the $\rho_{e}$ reconstructions. Alternatively, by comparing Fig. 1 (g,h) with Fig. 1 (a,b), we see that SIRZ-3 reduces the bias in $Z_{e}$ and $\rho_{e}$ estimates.

In Fig. 2 and 3, we compare the performance of SIRZ-2 and SIRZ-3 using the relative error metric (equation (14)). We compare the relative errors in $Z_{e}$ and $\rho_{e}$ for various materials with $Z_{e}$ in the range of $5-30$ and maximum low-energy attenuation values, denoted as $y_{L}^{max}$, in the range of $0.5-5.5$. Note that $y_{L}^{max}$ is the maximum value of low-energy attenuation (modeled using (6)) in the absence of noise. The rows of the heat map are sorted according to ascending values of $Z_{e}$. In Fig. 2 and 3, pure element materials are specified by their chemical formula, PVC ($Z_{e}^{GT}=14.07$)³³3$Z^{GT}_{e}$ denotes the ground-truth value for $Z_{e}$. denotes polyvinyl chloride, RbBr is $19\%$ RbBr salt solution in water ($Z_{e}^{GT}=20.05,\rho^{GT}_{e}=0.614\times 10^{-3}$electrons$\times$mol/mm³), Ti64 ($Z_{e}^{GT}=21.67$) is an alloy of Al ($2.8\%$), Ti ($23.9\%$), and V, Ti5553 ($Z_{e}^{GT}=23.73$) is an alloy of Al ($3.6\%$), Ti ($32.8\%$), V ($1.9\%$), Cr ($1.1\%$), and Mo, and brass is an alloy of Cu ($1.5\%$) and Zn.

The reduction in $Z_{e}$ relative error obtained using SIRZ-3 is quantified in Fig. 2 (c), which shows the difference between the absolute values of the relative errors using SIRZ-2 (Fig. 2 (a)) and SIRZ-3 (Fig. 2 (b)). A positive value in Fig. 2 (c) indicates that the SIRZ-2 relative error is larger in magnitude than the SIRZ-3 relative error. For LDPE, we observe that the SIRZ-3 absolute⁴⁴4Henceforth, relative error refers to the absolute value of the relative error. relative error is more than $1.5\%$ lower than the SIRZ-2 relative error when $y_{L}^{max}\leq 2.5$. Beyond LDPE, at low values for $Z_{e}$ of up to $Z_{e}=20$ (calcium), we see that both SIRZ-2 and SIRZ-3 errors are within $1\%$. Hence, there does not seem to be a benefit to using SIRZ-3 over SIRZ-2 since all errors are below $1\%$. With RbBr salt solution, we also do not see any apparent advantage of SIRZ-3 over SIRZ-2. For titanium ($Z_{e}=22$) and its alloys, we observe a marginal improvement in $Z_{e}$ performance with SIRZ-3 ($\approx 0\%-2\%$ reduction in relative error). For materials such as iron (Fe), copper (Cu), brass, and zinc (Zn) with $Z_{e}\geq 26$, the SIRZ-2 relative errors rise beyond $14\%$ when $y_{L}^{max}\geq 3.5$. However, we observe that the SIRZ-3 relative error is always less than $2.5\%$ for these materials. From Fig. 2, we observe a trend where the performance of SIRZ-3 over SIRZ-2 improves with increasing values for both $y_{L}^{max}$ and $Z_{e}$ when $Z_{e}$ is greater than $\approx 22$. The $Z_{e}$ relative errors with SIRZ-3 are almost always within $2.5\%$ (except RbBr) irrespective of the $Z_{e}$ value. However, the SIRZ-2 relative error increases to more than $25\%$ as the value for $Z_{e}$ and $y_{L}^{max}$ is increased. In Fig. 2, SIRZ-3 improves upon SIRZ-2 in $80\%$ of the simulated cases where the absolute value of the $Z_{e}$ relative errors with either SIRZ-2 or SIRZ-3 is greater than $1\%$, .

The reduction in $\rho_{e}$ relative error obtained using SIRZ-3 is quantified in Fig. 3 (c), which shows the difference between the absolute relative errors using SIRZ-2 (Fig. 3 (a)) and SIRZ-3 (Fig. 3 (b)). A positive value in Fig. 3 (c) indicates that the SIRZ-2 relative error is larger in magnitude than the SIRZ-3 relative error. At low values for $Z_{e}$ of up to PVC, there does not seem to be a benefit to using SIRZ-3 since the relative errors in $\rho_{e}$ are always under $1\%$. For calcium, SIRZ-3 reduces the $\rho_{e}$ relative error by more than $1\%$ when $y^{max}_{L}\geq 2.5$. SIRZ-3 does not seem to offer any improvements for RbBr salt solution. For Ti and its alloys, the reduction in SIRZ-3 relative error when compared to SIRZ-2 is in the range of $0.5-4\%$ when $y_{L}^{max}\geq 2.5$. For materials such as Fe, Cu, Brass, and Zn with $Z_{e}\geq 26$, we observe more than $20\%$ reduction in relative error using SIRZ-3 when $y_{L}^{max}\geq 2.5$. From Fig. 3, we observe a general trend where the $\rho_{e}$ performance of SIRZ-3 over SIRZ-2 improves with increasing values for $y_{L}^{max}$ and $Z_{e}$ when $Z_{e}$ is greater than $\approx 22$. While the SIRZ-2 relative errors in $\rho_{e}$ reach very large values of more than $40\%$, the SIRZ-3 relative errors in $\rho_{e}$ are almost always under $5\%$ (except RbBr). In Fig. 3, SIRZ-3 improves upon SIRZ-2 in $87\%$ of the simulated cases where the absolute value of the $\rho_{e}$ relative errors with either SIRZ-2 or SIRZ-3 is greater than $1\%$.

We presented the SIRZ-3 method for reconstruction of the effective atomic number ($Z_{e}$) and electron density ($\rho_{e}$) of an object imaged using X-ray dual-energy computed tomography (DECT). We presented a differentiable forward model that expresses the DECT data as a direct analytical function of the unknown $\left(Z_{e},\rho_{e}\right)$. SIRZ-3 uses numerical optimization to solve for $\left(Z_{e},\rho_{e}\right)$ such that the output of the forward model is close to the measured data. In general, compared to SIRZ-2, the magnitude of performance improvement with SIRZ-3 increases as a function of increasing $Z_{e}$. For materials such as iron, copper, brass, and zinc with $Z_{e}\geq 26$, we show that SIRZ-3 provides more than $10\%$ reduction in $Z_{e}$ relative error and more than $20\%$ reduction in $\rho_{e}$ relative error when the maximum low-energy attenuation is greater than $2.5$. The large errors for SIRZ-2 at high $Z_{e}$ is due to approximations made in the dual-energy decomposition and $\rho_{e}/Z_{e}$ conversion steps. Due to the absence of these approximations, SIRZ-3 performs better than SIRZ-2. We tested materials with $Z_{e}$ values in the range of $5-30$ and low-energy attenuation values in the range of $0.5-5.5$. During $Z_{e}$ reconstruction, SIRZ-3 improves upon SIRZ-2 in $80\%$ of our simulated cases where the absolute value of the $Z_{e}$ relative errors with either SIRZ-2 or SIRZ-3 is greater than $1\%$. During $\rho_{e}$ reconstruction, SIRZ-3 improves upon SIRZ-2 in $87\%$ of our simulated cases where the absolute value of the $\rho_{e}$ relative errors with either SIRZ-2 or SIRZ-3 is greater than $1\%$.

LLNL-CONF-832227. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. LDRD 21-FS-013 was used for this project.