DeepAI
Log In Sign Up

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

03/31/2022
by   K. Aditya Mohan, et al.
Lawrence Livermore National Laboratory
Arizona State University
0

For material identification, characterization, and quantification, it is useful to estimate system-independent material properties that do not depend on the detailed specifications of the X-ray computed tomography (CT) system such as spectral response. System independent rho-e and Z-e (SIRZ) refers to a suite of methods for estimating the system independent material properties of electron density, rho-e, and effective atomic number, Z-e, of an object scanned using dual-energy X-ray CT (DECT). The current state-of-the-art approach, SIRZ-2, makes certain approximations that lead to inaccurate estimates for large atomic numbered (Z-e) materials. In this paper, we present an extension, SIRZ-3, which iteratively reconstructs the unknown rho-e and Z-e while avoiding the limiting approximations made by SIRZ-2. Unlike SIRZ-2, this allows SIRZ-3 to accurately reconstruct rho-e and Z-e even at large values of Z-e. SIRZ-3 relies on the use of a new non-linear differentiable forward measurement model that expresses the DECT measurement data as a direct analytical function of rho-e and Z-e. Leveraging this new forward model, we use an iterative optimization algorithm to reconstruct (or solve for) rho-e and Z-e directly from the DECT data. Compared to SIRZ-2, we show that the magnitude of performance improvement using SIRZ-3 increases with increasing values for Z-e.

READ FULL TEXT VIEW PDF

page 5

page 6

page 7

01/06/2020

Uniqueness criteria in multi-energy CT

Multi-Energy Computed Tomography (ME-CT) is a medical imaging modality a...
05/02/2019

A Splitting-Based Iterative Algorithm for GPU-Accelerated Statistical Dual-Energy X-Ray CT Reconstruction

When dealing with material classification in baggage at airports, Dual-E...
09/29/2021

An inversion algorithm for P-functions with applications to Multi-energy CT

Multi-energy computed tomography (ME-CT) is an x-ray transmission imagin...
12/21/2021

ADJUST: A Dictionary-Based Joint Reconstruction and Unmixing Method for Spectral Tomography

Advances in multi-spectral detectors are causing a paradigm shift in X-r...
03/25/2021

Regularization by Denoising Sub-sampled Newton Method for Spectral CT Multi-Material Decomposition

Spectral Computed Tomography (CT) is an emerging technology that enables...
04/21/2022

Parametric Level-sets Enhanced To Improve Reconstruction (PaLEnTIR)

In this paper, we consider the restoration and reconstruction of piecewi...

1 Introduction

X-ray 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 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 [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 SIRZ-2 approach that improved upon the original SIRZ to produce more accurate estimates for and . SIRZ-2 can be broadly split into four different steps:

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

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

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

  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 and 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 and .

  • Formulation of an error function that quantifies the deviation of measured data from the output of the forward model given the unknowns, .

  • Method for reconstructing from the measured data by minimizing the error function using the limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) solver [5, 7].

2 Background

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 , on the electron density, , and atomic number, , is given by,

(1)

where is the X-ray energy, is the total X-ray electronic cross-section.[2] Equation (1) is only defined for integer valued that correspond to pure elements.

In order to represent compounds or mixtures of elements, the cross-section for non-integer 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 X-ray cross-section with superior system invariance properties. The for a compound or mixture is defined as,

(3)

where is the tabulated[4] electron cross-section 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 X-ray energy range of interest (30-200keV). 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 X-ray energy at the energy bin. The linear projections (line-integration) 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 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, .

For DECT, the complete forward model is,

(6)
(7)

where is the spectral density for X-rays at the low-energy spectrum, is the spectral density at the high-energy spectrum, and and 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 with respect to . This derivative is given by,

(8)

where and is the floor function.

4 Sirz-3: Reconstruction by Optimization

SIRZ-3 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 . SIRZ-3 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 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,

(9)

where and are the attenuation space measurements at the low-energy and high-energy 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 memory-efficient sparse representation that only stores the non-zero values and its locations. We used LTT [3] to generate . Then, the reconstruction of SIRZ-3 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 z-score normalized versions of

and . Let and

be the mean and standard deviation of the SIRZ-2 reconstructed voxel values for

. Similarly, let and be the mean and standard deviation of the SIRZ-2 reconstructed voxel values for . Then, we solve,

(12)

where we use vector-scalar multiplication111Vector-scalar multiplication is point wise multiplication of a scalar to every element of a vector. and vector-scalar addition222Vector-scalar addition is point wise addition of a scalar to every element of a vector. within the error function . The SIRZ-3 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 L-BFGS [7, 5] optimization algorithm. Upon convergence and using (13), L-BFGS 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 SIRZ-2 estimates for and to initialize the optimization in (12).

Ginsimsz[true]Ginheight=1.15in, keepaspectratio=true

(a) Ground-Truth (b) Ground-Truth (c) Low-Energy Sinogram (d) High-Energy Sinogram
(e) SIRZ-2 (f) SIRZ-2 (g) SIRZ-3 (h) SIRZ-3
Relative Error: Relative Error: Relative Error: Relative Error:
Relative RMSE: Relative RMSE: Relative RMSE: Relative RMSE:
Figure 1: () reconstruction of a mm diameter copper (Cu) disc with a maximum attenuation of at the low-energy X-ray spectrum. (a, b) show the ground-truth images. (c, d) show the sinograms at the low-energy and high-energy spectra. The axis in (c,d) is along the view dimension. (e, f) and (g, h) show the SIRZ-2 and SIRZ-3 reconstructions respectively. Compared to SIRZ-2, SIRZ-3 reduces both the relative error (equation (14)) and relative RMSE (equation (15)) in and reconstructions. The high relative error associated with SIRZ-2 is also observed by visually comparing (e, f) with (a,b). The noisy ring-like artifacts around the reconstruction of the copper disk is due to the near zero density of the surrounding air. The unit for is electronsmol/mm.

Comparison of relative errors (%) in as defined in equation (14)

(a) SIRZ-2 Relative Error, (b) SIRZ-3 Relative Error, (c) Difference Relative Error,
Figure 2: Heat map of the relative error (equation (14)) in reconstruction of for various materials and maximum low-energy attenuation values, . (a) and (b) show the relative errors for SIRZ-2 and SIRZ-3 respectively. (c) shows the difference between the absolute values of the SIRZ-2 relative error in (a) and SIRZ-3 relative error in (b). The rows in the heat map are sorted according to the value. In general, compared to SIRZ-2, the reduction in relative error of using SIRZ-3 improves with increasing values for . When the absolute value of the relative errors with either SIRZ-2 or SIRZ-3 is greater than , SIRZ-3 improves upon SIRZ-2 in of the simulated cases in the heat map.

Comparison of relative errors (%) in as defined in equation (14)

(a) SIRZ-2 Relative Error, (b) SIRZ-3 Relative Error, (c) Difference Relative Error,
Figure 3: Heat map of the relative error (equation (14)) in reconstruction of for various materials and maximum low-energy attenuation values. (a) and (b) show the relative errors (equation (14)) for SIRZ-2 and SIRZ-3 respectively. (c) shows the difference between the absolute values of the SIRZ-2 relative error in (a) and SIRZ-3 relative error in (b). The rows in the heat map are sorted according to the value. In general, compared to SIRZ-2, the reduction in relative error of using SIRZ-3 improves with increasing values for . When the absolute value of the relative errors with either SIRZ-2 or SIRZ-3 is greater than , SIRZ-3 improves upon SIRZ-2 in of the simulated cases in the heat map.

5 Results

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, , in the range of 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 . 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 and object thicknesses.

We simulated fan-beam X-ray DECT data at different view angles and pixel bins. The low-energy X-ray spectrum is generated using a kV bremsstrahlung X-ray source and a mm Aluminum filter. The high-energy X-ray spectrum is generated using a kV bremsstrahlung X-ray source, a mm aluminum filter, and a mm copper filter. DECT data is simulated for a circular object that occupies a cross-sectional pixel area within the full field-of-view of . The source-to-object (SOD) distance and source-to-detector (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 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 (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 [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 SIRZ-2 or SIRZ-3 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 ground-truth 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 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 and a maximum attenuation of at the low-energy spectrum. Fig. 1 (a, b) show the ground-truth 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 and under estimates the reconstructions. Alternatively, by comparing Fig. 1 (g,h) with Fig. 1 (a,b), we see that SIRZ-3 reduces the bias in and 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 and for various materials with in the range of and maximum low-energy attenuation values, denoted as , in the range of . Note that 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 . In Fig. 2 and 3, pure element materials are specified by their chemical formula, PVC ()333 denotes the ground-truth 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 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 absolute444Henceforth, relative error refers to the absolute value of the relative error. relative error is more than lower than the SIRZ-2 relative error when . Beyond LDPE, at low values for of up to (calcium), we see that both SIRZ-2 and SIRZ-3 errors are within . Hence, there does not seem to be a benefit to using SIRZ-3 over SIRZ-2 since all errors are below . With RbBr salt solution, we also do not see any apparent advantage of SIRZ-3 over SIRZ-2. For titanium () and its alloys, we observe a marginal improvement in performance with SIRZ-3 ( reduction in relative error). For materials such as iron (Fe), copper (Cu), brass, and zinc (Zn) with , the SIRZ-2 relative errors rise beyond when . However, we observe that the SIRZ-3 relative error is always less than 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 and when is greater than . The relative errors with SIRZ-3 are almost always within (except RbBr) irrespective of the value. However, the SIRZ-2 relative error increases to more than as the value for and is increased. In Fig. 2, SIRZ-3 improves upon SIRZ-2 in of the simulated cases where the absolute value of the relative errors with either SIRZ-2 or SIRZ-3 is greater than , .

The reduction in 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 of up to PVC, there does not seem to be a benefit to using SIRZ-3 since the relative errors in are always under . For calcium, SIRZ-3 reduces the relative error by more than when . 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 when . For materials such as Fe, Cu, Brass, and Zn with , we observe more than reduction in relative error using SIRZ-3 when . From Fig. 3, we observe a general trend where the performance of SIRZ-3 over SIRZ-2 improves with increasing values for and when is greater than . While the SIRZ-2 relative errors in reach very large values of more than , the SIRZ-3 relative errors in are almost always under (except RbBr). In Fig. 3, SIRZ-3 improves upon SIRZ-2 in of the simulated cases where the absolute value of the relative errors with either SIRZ-2 or SIRZ-3 is greater than .

6 Conclusions

We presented the SIRZ-3 method for reconstruction of the effective atomic number () and electron density () 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 . SIRZ-3 uses numerical optimization to solve for 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 . For materials such as iron, copper, brass, and zinc with , we show that SIRZ-3 provides more than reduction in relative error and more than reduction in relative error when the maximum low-energy attenuation is greater than . The large errors for SIRZ-2 at high is due to approximations made in the dual-energy decomposition and conversion steps. Due to the absence of these approximations, SIRZ-3 performs better than SIRZ-2. We tested materials with values in the range of and low-energy attenuation values in the range of . During reconstruction, SIRZ-3 improves upon SIRZ-2 in of our simulated cases where the absolute value of the relative errors with either SIRZ-2 or SIRZ-3 is greater than . During reconstruction, SIRZ-3 improves upon SIRZ-2 in of our simulated cases where the absolute value of the relative errors with either SIRZ-2 or SIRZ-3 is greater than .

7 Acknowledgments

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.

References

  • [1] S. G. Azevedo, H. E. Martz, M. B. Aufderheide, W. D. Brown, K. M. Champley, J. S. Kallman, G. P. Roberson, D. Schneberk, I. M. Seetho, and J. A. Smith (2016) System-independent characterization of materials using dual-energy computed tomography. IEEE Transactions on Nuclear Science 63 (1), pp. 341–350. External Links: Document Cited by: §1, §2.
  • [2] K. M. Champley, S. G. Azevedo, I. M. Seetho, S. M. Glenn, L. D. McMichael, J. A. Smith, J. S. Kallman, W. D. Brown, and H. E. Martz (2019) Method to extract system-independent material properties from dual-energy x-ray CT. IEEE Transactions on Nuclear Science 66 (3), pp. 674–686. External Links: Document Cited by: §1, §1, §2, §2.
  • [3] K. M. Champley, T. M. Willey, H. Kim, K. Bond, S. M. Glenn, J. A. Smith, J. S. Kallman, W. D. Brown, I. M. Seetho, L. Keene, S. G. Azevedo, L. D. McMichael, G. Overturf, and H. E. Martz (2022) Livermore tomography tools: Accurate, fast, and flexible software for tomographic science. NDT & E International 126, pp. 102595. External Links: ISSN 0963-8695, Link, Document Cited by: §2, §4, §5.
  • [4] D. E. Cullen, J. H. Hubbell, and L. Kissel (1997) EPDL97: the evaluated photo data library 97 version. Lawrence Livermore National Lab., CA (United States). External Links: Document Cited by: §2.
  • [5] D. C. Liu and J. Nocedal (1989) On the limited memory BFGS method for large scale optimization. Mathematical programming 45 (1), pp. 503–528. Cited by: 3rd item, §4.
  • [6] H. E. Martz, C. M. Logan, D. J. Schneberk, and P. J. Shull (2016) X-ray imaging: fundamentals, industrial techniques and applications. CRC Press. Cited by: §1.
  • [7] J. Nocedal (1980) Updating quasi-newton matrices with limited storage. Mathematics of computation 35 (151), pp. 773–782. Cited by: 3rd item, §4.
  • [8] H. M. Shi and D. Mudigere PyTorch-LBFGS: a PyTorch implementation of L-BFGS. https://github.com/hjmshi/PyTorch-LBFGS. Cited by: §4.