1 Introduction
Xray computed tomography (CT) is widely used to reconstruct the spatially varying linear attenuation coefficient (LAC) of an object [6]. The spectral response function that is dependent on the Xray source and the detector is used to quantify the variation in Xray intensity over the broad range of emitted Xray energies. Most CT systems have differing spectral response functions that causes the reconstructed LAC to be dependent on the spectral response of the Xray scanner. Hence, it is not possible to compare the LAC values reconstructed using different Xray scanners. Also, the estimated LAC does not provide direct information on the density and atomic composition of an object [2].
Dualenergy Xray 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 [1] to estimate the system independent properties of electron density, , and effective atomic number, . The electron density, , is a direct measure of the object’s density and the effective atomic number, , is useful to predict the atomic composition. In Champley at al.[2], the authors presented the SIRZ2 approach that improved upon the original SIRZ to produce more accurate estimates for and . SIRZ2 can be broadly split into four different steps:

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

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

Filtered backprojection to reconstruct the LAC values at the two SMB energies,

Estimation of electron density and effective atomic number from the LACs at the two SMB energies.
The second and fourth steps of SIRZ2 are based on simplifying approximations that limit its accuracy.
In this paper, we present a new approach called SIRZ3 to reconstruct the system independent properties of and while avoiding the approximations made in SIRZ2. SIRZ3 replaces the steps 2, 3, and 4 of SIRZ2 using a single iterative optimization algorithm. Features of the SIRZ3 method include:

Nonlinear differentiable forward model that expresses the measured DECT data as a direct function of the object’s spatial distribution of and .

Formulation of an error function that quantifies the deviation of measured data from the output of the forward model given the unknowns, .
2 Background
The LAC quantifies the magnitude of Xray attenuation within an object. The LAC varies as a function of the position inside the object, density, atomic composition, and Xray energy. The dependence of LAC, denoted as , on the electron density, , and atomic number, , is given by,
(1) 
where is the Xray energy, is the total Xray electronic crosssection.[2] Equation (1) is only defined for integer valued that correspond to pure elements.
In order to represent compounds or mixtures of elements, the crosssection for noninteger is defined as follows [1],
(2) 
where , and . Note that for a given compound or mixture is not a fundamental physical constant and there are several competing techniques to define . We adopt the definition of from Champley et al. [2] since it is a function of the Xray crosssection with superior system invariance properties. The for a compound or mixture is defined as,
(3) 
where is the tabulated[4] electron crosssection of the material, is the spectral response at energy , and is the areal electron density. This definition in (3) is demonstrated to be relatively insensitive to and the Xray energy range of interest (30200keV). For additional details on the definition of and its estimation, the reader may refer to Champley et al. [3].
3 Forward Model
In this section, we formulate a forward model that mathematically expresses the measured DECT data as a function of the object properties, and . This forward model is necessary to solve the inverse problem of reconstructing the object’s and .
The forward model is derived in discrete coordinate space. Let and represent the and values respectively at the voxel inside the object. Then, the LAC at voxel is given by , where from (2) and is the Xray energy at the energy bin. The linear projections (lineintegration) of the LAC is expressed as,
(4) 
where is the element along the row and column of the matrix . At the detector, the measurement is the average transmission function weighted by the spectral density function . Hence, the effective transmission is given by,
(5) 
where is the spectral density at the energy bin ( is a discretized version of used in (3)) such that . The spectral density function quantifies both the spectral density of the Xray 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, .
For DECT, the complete forward model is,
(6)  
(7) 
where is the spectral density for Xrays at the lowenergy spectrum, is the spectral density at the highenergy spectrum, and and are the forward model predictions for the lowenergy and highenergy measurements respectively. To ensure differentiability of the forward model in (6) and (7), we need analytical derivatives for with respect to . This derivative is given by,
(8) 
where and is the floor function.
4 Sirz3: Reconstruction by Optimization
SIRZ3 is a framework to solve the problem of reconstructing and from measurements and using a forward model that accurately defines an analytical relation between and for all and . SIRZ3 adopts an iterative approach, where the solution for is iteratively refined such that the output of the forward model (whose input is the current estimate for and ) gets progressively closer to the measurements . In addition to a forward model, the solution for 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 SIRZ3 is dependent on the accuracy of this analytical relation that mathematically describes the measurement process of our Xray system. Future efforts to improve SIRZ3 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,
(9) 
where and are the attenuation space measurements at the lowenergy and highenergy spectra respectively, and and are the forward model predictions from (6) and (7) respectively. The weights for the penalty terms are given by and , where is the number of sinogram pixels summed over all the views. Substituting (6) and (7) in (9), we get,
(10) 
where
is a vector of all
voxel values, is a vector of all voxel values. Note that the projection matrix entries are stored in a memoryefficient sparse representation that only stores the nonzero values and its locations. We used LTT [3] to generate . Then, the reconstruction of SIRZ3 is obtained by solving,(11) 
where and are reconstructions of and respectively, and are the lower and upper limits for respectively, and and are the lower and upper limits for respectively. In this paper, we set , , , and electronsmol/mm.
In order to improve conditioning and ensure fast convergence, we optimize zscore normalized versions of
and . Let andbe the mean and standard deviation of the SIRZ2 reconstructed voxel values for
. Similarly, let and be the mean and standard deviation of the SIRZ2 reconstructed voxel values for . Then, we solve,(12) 
where we use vectorscalar multiplication^{1}^{1}1Vectorscalar multiplication is point wise multiplication of a scalar to every element of a vector. and vectorscalar addition^{2}^{2}2Vectorscalar addition is point wise addition of a scalar to every element of a vector. within the error function . The SIRZ3 reconstruction is given by,
(13) 
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 with respect to and . We solve equation (12) using the LBFGS [7, 5] optimization algorithm. Upon convergence and using (13), LBFGS yields the reconstruction of the electron density, , and the effective atomic number, . We use the LBFGS implementation by Shi et al. [8]. For LBFGS, we use a history size of , weak Wolfe line search, and a convergence criteria that stops the LBFGS iterations when the percentage change in both and goes below for 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 SIRZ2 estimates for and to initialize the optimization in (12).Ginsimsz[true]Ginheight=1.15in, keepaspectratio=true
(a) GroundTruth  (b) GroundTruth  (c) LowEnergy Sinogram  (d) HighEnergy Sinogram 
(e) SIRZ2  (f) SIRZ2  (g) SIRZ3  (h) SIRZ3 
Relative Error:  Relative Error:  Relative Error:  Relative Error: 
Relative RMSE:  Relative RMSE:  Relative RMSE:  Relative RMSE: 
Comparison of relative errors (%) in as defined in equation (14)
(a) SIRZ2 Relative Error,  (b) SIRZ3 Relative Error,  (c) Difference Relative Error, 

Comparison of relative errors (%) in as defined in equation (14)
(a) SIRZ2 Relative Error,  (b) SIRZ3 Relative Error,  (c) Difference Relative Error, 
5 Results
Using simulated data, we compared the performance of SIRZ2 and SIRZ3 for a wide range of materials and object thicknesses. We chose various materials including pure elements, compounds, and mixtures with effective atomic numbers, , in the range of since they are among the most commonly occurring materials that are sensitive to Xrays in our chosen energy range. Our choice of object’s thicknesses was such that the maximum attenuation value for the lowenergy DECT spectrum was in the range of . Outside this attenuation range, it is typically difficult to extract information from Xray 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 and object thicknesses.
We simulated fanbeam Xray DECT data at different view angles and pixel bins. The lowenergy Xray spectrum is generated using a kV bremsstrahlung Xray source and a mm Aluminum filter. The highenergy Xray spectrum is generated using a kV bremsstrahlung Xray source, a mm aluminum filter, and a mm copper filter. DECT data is simulated for a circular object that occupies a crosssectional pixel area within the full fieldofview of . The sourcetoobject (SOD) distance and sourcetodetector (SDD) distance were chosen to be mm and mm respectively. The diameter of the object was adjusted such that we obtained maximum attenuation values of to in increments of for the lowenergy spectrum. For each attenuation value, given the atomic composition of the object’s material, we computed the diameter that produces the desired lowenergy attenuation along the Xray path passing through the center of the object. Then, the pixel size for the object was obtained by dividing this diameter by (number of pixels occupied by the object along any direction). We assume a PerkinElmer (PE) detector for the detector spectral response. For generating simulated data, we used analytic raytracing through geometric solids that is implemented in LTT [3]. We add noise to the measurements in transmission space.
For quantitative comparison of and 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 or as,
(14) 
where is used to denote either SIRZ2 or SIRZ3 methods, denotes either or , is the set of all pixel indices in the interior of the object that excludes the object boundaries, is the cardinality (number of elements) of the set , and is the groundtruth value for . The second error metric is the relative root mean squared error (RMSE), which is defined as,
(15) 
The relative RMSE metric combines both the bias and standard deviation of the difference between the reconstruction and the groundtruth.
In Fig. 1, we compare the performance of SIRZ2 and SIRZ3 for a disc of copper (Cu) with a diameter of and a maximum attenuation of at the lowenergy spectrum. Fig. 1 (a, b) show the groundtruth images and Fig. 1 (c, d) show the lowenergy and highenergy sinograms. SIRZ2 and SIRZ3 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 SIRZ3 when compared to SIRZ2. By comparing Fig. 1 (e,f) with Fig. 1 (a,b), we see that SIRZ2 over estimates the and under estimates the reconstructions. Alternatively, by comparing Fig. 1 (g,h) with Fig. 1 (a,b), we see that SIRZ3 reduces the bias in and estimates.
In Fig. 2 and 3, we compare the performance of SIRZ2 and SIRZ3 using the relative error metric (equation (14)). We compare the relative errors in and for various materials with in the range of and maximum lowenergy attenuation values, denoted as , in the range of . Note that is the maximum value of lowenergy attenuation (modeled using (6)) in the absence of noise. The rows of the heat map are sorted according to ascending values of . In Fig. 2 and 3, pure element materials are specified by their chemical formula, PVC ()^{3}^{3}3 denotes the groundtruth value for . denotes polyvinyl chloride, RbBr is RbBr salt solution in water (electronsmol/mm), Ti64 () is an alloy of Al (), Ti (), and V, Ti5553 () is an alloy of Al (), Ti (), V (), Cr (), and Mo, and brass is an alloy of Cu () and Zn.
The reduction in relative error obtained using SIRZ3 is quantified in Fig. 2 (c), which shows the difference between the absolute values of the relative errors using SIRZ2 (Fig. 2 (a)) and SIRZ3 (Fig. 2 (b)). A positive value in Fig. 2 (c) indicates that the SIRZ2 relative error is larger in magnitude than the SIRZ3 relative error. For LDPE, we observe that the SIRZ3 absolute^{4}^{4}4Henceforth, relative error refers to the absolute value of the relative error. relative error is more than lower than the SIRZ2 relative error when . Beyond LDPE, at low values for of up to (calcium), we see that both SIRZ2 and SIRZ3 errors are within . Hence, there does not seem to be a benefit to using SIRZ3 over SIRZ2 since all errors are below . With RbBr salt solution, we also do not see any apparent advantage of SIRZ3 over SIRZ2. For titanium () and its alloys, we observe a marginal improvement in performance with SIRZ3 ( reduction in relative error). For materials such as iron (Fe), copper (Cu), brass, and zinc (Zn) with , the SIRZ2 relative errors rise beyond when . However, we observe that the SIRZ3 relative error is always less than for these materials. From Fig. 2, we observe a trend where the performance of SIRZ3 over SIRZ2 improves with increasing values for both and when is greater than . The relative errors with SIRZ3 are almost always within (except RbBr) irrespective of the value. However, the SIRZ2 relative error increases to more than as the value for and is increased. In Fig. 2, SIRZ3 improves upon SIRZ2 in of the simulated cases where the absolute value of the relative errors with either SIRZ2 or SIRZ3 is greater than , .
The reduction in relative error obtained using SIRZ3 is quantified in Fig. 3 (c), which shows the difference between the absolute relative errors using SIRZ2 (Fig. 3 (a)) and SIRZ3 (Fig. 3 (b)). A positive value in Fig. 3 (c) indicates that the SIRZ2 relative error is larger in magnitude than the SIRZ3 relative error. At low values for of up to PVC, there does not seem to be a benefit to using SIRZ3 since the relative errors in are always under . For calcium, SIRZ3 reduces the relative error by more than when . SIRZ3 does not seem to offer any improvements for RbBr salt solution. For Ti and its alloys, the reduction in SIRZ3 relative error when compared to SIRZ2 is in the range of when . For materials such as Fe, Cu, Brass, and Zn with , we observe more than reduction in relative error using SIRZ3 when . From Fig. 3, we observe a general trend where the performance of SIRZ3 over SIRZ2 improves with increasing values for and when is greater than . While the SIRZ2 relative errors in reach very large values of more than , the SIRZ3 relative errors in are almost always under (except RbBr). In Fig. 3, SIRZ3 improves upon SIRZ2 in of the simulated cases where the absolute value of the relative errors with either SIRZ2 or SIRZ3 is greater than .
6 Conclusions
We presented the SIRZ3 method for reconstruction of the effective atomic number () and electron density () of an object imaged using Xray dualenergy computed tomography (DECT). We presented a differentiable forward model that expresses the DECT data as a direct analytical function of the unknown . SIRZ3 uses numerical optimization to solve for such that the output of the forward model is close to the measured data. In general, compared to SIRZ2, the magnitude of performance improvement with SIRZ3 increases as a function of increasing . For materials such as iron, copper, brass, and zinc with , we show that SIRZ3 provides more than reduction in relative error and more than reduction in relative error when the maximum lowenergy attenuation is greater than . The large errors for SIRZ2 at high is due to approximations made in the dualenergy decomposition and conversion steps. Due to the absence of these approximations, SIRZ3 performs better than SIRZ2. We tested materials with values in the range of and lowenergy attenuation values in the range of . During reconstruction, SIRZ3 improves upon SIRZ2 in of our simulated cases where the absolute value of the relative errors with either SIRZ2 or SIRZ3 is greater than . During reconstruction, SIRZ3 improves upon SIRZ2 in of our simulated cases where the absolute value of the relative errors with either SIRZ2 or SIRZ3 is greater than .
7 Acknowledgments
LLNLCONF832227. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DEAC5207NA27344. LDRD 21FS013 was used for this project.
References
 [1] (2016) Systemindependent characterization of materials using dualenergy computed tomography. IEEE Transactions on Nuclear Science 63 (1), pp. 341–350. External Links: Document Cited by: §1, §2.
 [2] (2019) Method to extract systemindependent material properties from dualenergy xray CT. IEEE Transactions on Nuclear Science 66 (3), pp. 674–686. External Links: Document Cited by: §1, §1, §2, §2.
 [3] (2022) Livermore tomography tools: Accurate, fast, and flexible software for tomographic science. NDT & E International 126, pp. 102595. External Links: ISSN 09638695, Link, Document Cited by: §2, §4, §5.
 [4] (1997) EPDL97: the evaluated photo data library 97 version. Lawrence Livermore National Lab., CA (United States). External Links: Document Cited by: §2.
 [5] (1989) On the limited memory BFGS method for large scale optimization. Mathematical programming 45 (1), pp. 503–528. Cited by: 3rd item, §4.
 [6] (2016) Xray imaging: fundamentals, industrial techniques and applications. CRC Press. Cited by: §1.
 [7] (1980) Updating quasinewton matrices with limited storage. Mathematics of computation 35 (151), pp. 773–782. Cited by: 3rd item, §4.
 [8] PyTorchLBFGS: a PyTorch implementation of LBFGS. https://github.com/hjmshi/PyTorchLBFGS. Cited by: §4.