X-ray Spectral Estimation using Dictionary Learning

TL;DR

This paper introduces a dictionary learning-based spectral estimation method for X-ray CT systems that accurately estimates spectral response without initial guesses, improving over existing methods.

Contribution

The paper presents a novel dictionary-based spectral estimation approach using MAP estimation with sparsity and simplex constraints, outperforming state-of-the-art methods.

Findings

Accurate spectral estimates achieved without initial guesses.

Method outperforms existing techniques in cross-validation experiments.

Demonstrated effectiveness on real ALS datasets.

Abstract

As computational tools for X-ray computed tomography (CT) become more quantitatively accurate, knowledge of the source-detector spectral response is critical for quantitative system-independent reconstruction and material characterization capabilities. Directly measuring the spectral response of a CT system is hard, which motivates spectral estimation using transmission data obtained from a collection of known homogeneous objects. However, the associated inverse problem is ill-conditioned, making accurate estimation of the spectrum challenging, particularly in the absence of a close initial guess. In this paper, we describe a dictionary-based spectral estimation method that yields accurate results without the need for any initial estimate of the spectral response. Our method utilizes a MAP estimation framework that combines a physics-based forward model along with an $L_0$ sparsity constraint and a simplex constraint on the dictionary coefficients. Our method uses a greedy support selection method and a new pair-wise iterated coordinate descent method to compute the above estimate. We demonstrate that our dictionary-based method outperforms a state-of-the-art method as shown in a cross-validation experiment on four real datasets collected at beamline 8.3.2 of the Advanced Light Source (ALS).

Tables

Table 1. Table 1 : Setup for X-ray Scanning

Projection Geometry:	Parallel beam geometry
Source filter:	2 $mm$ Silicon
Scintillator:	50 $\mu m$ $L{u}_{3}A{l}_5{O}_{12}$
Max Energy:	100 KeV
Views Spanning:	Equi-spaced in $[0,2\pi ]$
Detector Pixel Size:	0.00065 mm
Nviews $\times$ Nrows $\times$ Ncolumns:	$2625\times 100\times 2560$
Sample-detector distance:	$0.3$ mm

Table 2. Table 2 : Dictionary Generation List

Filter response
Material	Thickness Range $mm$	Step $mm$	# of responses
$Al$	0.1$\sim$5.9	0.2	30
$Cu$	0.2$\sim$0.49	0.01	30
Detector response(Scintillator)
$L{u}_{3}A{l}_5{O}_{12}$	0.025$\sim$0.095	0.002	36

Table 3. Table 3 : Leave-One-Out Cross-Validation NRMSE

Case	Fit	Test	LSSE	DictSE
1	$Ti,Mg,Al$	$V$	0.0331	0.0315
2	$V,Mg,Al$	$Ti$	0.0624	0.0242
3	$V,Ti,Al$	$Mg$	0.0343	0.0122
4	$V,Ti,Mg$	$Al$	0.0483	0.0093

Equations

$$S(E)=S_{sr}(E)S_{ft}(E)S_{dt}(E),$$

$$I={\int }_0^{E_{\max }}S(E)\cdot \exp \left\{-{\int }_L\mu (E,r)dr\right\}dE,$$

$$y$$

$$y_i$$

$$Y=Fx+\tau ,$$

$$\text{s.t. }\Vert x{\Vert }_1=1,x_j\ge 0\forall j,$$

$$x=D\omega$$

$$Y=FD\omega +\tau .$$

$$\hat{\omega }=\arg \underset{\begin{array}{c}\omega \in \mathcal{S}\\ \Vert \omega {\Vert }_0\le C\end{array}}{\min }l(\omega ),$$

$$l^{\prime }(k,\beta |\hat{\omega })$$

$$=\frac{1}{2}{\Vert Y-\beta FD\hat{\omega }-\left(1-\beta \right)F{D}_{\ast ,k}\Vert }_{\Lambda }^2$$

$$\left(k^{\ast },{\beta }^{\ast }\right)$$

$${\beta }_k=\frac{e_k^T\Lambda (Y-F{D}_{\ast ,k})}{\Vert e_k{\Vert }_{\Lambda }^2}.$$

$${\alpha }^{\ast }\leftarrow \arg \underset{\alpha \in \left[-{\omega }_k,{\omega }_g\right]}{\min }\left\{l\left(\omega +\alpha ({ϵ}_k-{ϵ}_g)\right)\right\}.$$

$${\alpha }^{\ast }=\text{Clip}\left\{\frac{(Y-FD\omega {)}^{\mathrm{T}}\Lambda FD({ϵ}_k-{ϵ}_g)}{{\Vert FD({ϵ}_k-{ϵ}_g)\Vert }_{\Lambda }^2},[-{\omega }_k,{\omega }_g]\right\}$$

Full text

In this paper, we describe a dictionary-based spectral estimation method that yields accurate results without the need for any initial estimate of the spectral response. Our method utilizes a MAP estimation framework that combines a physics-based forward model along with an $L_{0}$ sparsity constraint and a simplex constraint on the dictionary coefficients. Our method uses a greedy support selection method and a new pair-wise iterated coordinate descent method to compute the above estimate. We demonstrate that our dictionary-based method outperforms a state-of-the-art method as shown in a cross-validation experiment on four real datasets collected at beamline 8.3.2 of the Advanced Light Source (ALS).††Document Release Number: LLNL-CONF-845171

**Index Terms— ** X-ray CT, spectral estimation, dictionary learning, inverse problem

1 Introduction

Non-destructive evaluation (NDE) is an increasingly important application of CT, which motivates improved capabilities for quantitative energy-independent reconstruction and precise material characterization. Such reconstruction techniques utilize physically accurate forward models that account for the polyenergetic nature of the X-ray source radiation and associated physical phenomena such as beam-hardening[1], rather than simpler models based on mono-energetic approximations. However, an accurate estimate of the CT system’s source-detector spectral response is a necessity for the above methods. For example, the method for reconstructing energy-independent material properties like effective atomic-number and electron density from dual-energy CT scans [2, 3, 4] and the method for tissue characterization in [5, 6] require a precise calibration of the source-detector spectral response.

Direct measurement of the spectral response of the whole X-Ray CT system is difficult since the detector response is hard to measure. Based on Beer–Lambert’s law[7], we can use a linear model for transmission measurements of objects with known dimensions and composition to reconstruct the discretized spectral response of a CT scanner, but this yields a highly ill-conditioned system. Champley et al. [4] use linear least-squares to do spectral estimation (LSSE) with constraints to enforce non-negativity. However, this method requires an accurate initial guess close to the true spectrum. Various regularization methods[8, 9] are used to overcome the issues introduced by the ill-conditioned nature of the problem. SVD-based algorithms have also been applied to the spectral estimation problem[10, 11, 12]. Sidky et al.[13] represent the spectrum as a linear combination of B-splines and use EM (Expectation–maximization) to find the solution. Zhao et al.[14] estimate the spectrum as a linear combination of six Monte Carlo model spectra. Liu et al.[15] introduce compressed sensing to estimate the spectrum. The central ideas of the above methods have a common theme: try to solve the ill-conditioned inverse problem either with regularized optimization or by introducing basis spectra to perform the estimation.

In this paper, we introduce a novel dictionary-based spectral estimation (DictSE) method that can efficiently reconstruct the overall spectral response of a CT system from transmission scans of multiple known objects, without the need for accurate initialization. We represent the unknown spectral response using an over-complete dictionary that accounts for vast combinations of different source spectra, filter attenuation characteristics, and detector energy-response models. We formulate the reconstruction problem as a MAP estimation framework that combines a linear beam-hardening forward model along with prior constraints. Specifically, we impose an $L_{0}$ sparsity constraint to limit the support for the spectrum representation and a simplex constraint to account for the bright-dark normalization of the transmission data. We present a novel iterative optimization strategy that alternates between support selection and pairwise iterative coordinate descent (ICD) update to find the optimal sparse representation of the spectrum. Finally, we demonstrate DictSE through a cross-validation experiment on four datasets collected at beamline 8.3.2 of ALS.

2 Reconstruction model

In this section, we describe the spectral estimation problem and our proposed solution in more detail. We use a physics-based model for X-ray transmission measurements and discretize the model in energy and by projection to produce a linear measurement model. Then we introduce a dictionary-based framework and express the reconstruction as a MAP estimation problem. Finally, we solve this dictionary-based MAP problem using support selection to enforce sparse coding and an ICD algorithm modified to enforce a simplex constraint on the coefficients of the selected dictionary elements.

2.1 Physics-based Model and Discretization

Fig. 1 illustrates the setup for X-ray spectral estimation. Using notation as in that figure, the spectral response of this CT system is the product of the X-ray source spectrum, filter response, and detector response, which yields the response function

[TABLE]

and our goal is to estimate $S(E)$.

We relate $S(E)$ to measurements by first defining $I$ and $I_{0}$ to be the measured intensity of the object scan (with sample) and the blank scan (without sample), respectively. Based on Beer–Lambert’s law, we have

[TABLE]

where $\mu(E,r)$ is linear attenuation coefficient (LAC) with units $mm^{-1}$. Likewise, $I_{0}$ is defined as in (2) with $\mu=0$. The (normalized) transmission through line $L$ is then defined as

[TABLE]

where $\bar{S}(E)=\frac{S(E)}{\int_{0}^{E_{\max}}S(E)dE}$ is the normalized response.

We assume each sample consists of solid rods made from known materials taken from a specified reference set $\Phi$. For each material $s\in\Phi$ and projection $i\in\{0,...,M-1\}$, we define $L_{i,s}$ to be the path length of the $i^{th}$ projection through the $s^{th}$ material and $\mu_{s}(E)$ to be the LAC of the $s$ material at energy $E$. Assuming that $I_{0}$ is independent of projection and noise-free, this yields the transmission for the $i^{th}$ projection as

[TABLE]

where $\tau_{i}$ is additive noise.

We discretize in energy by subdividing into non-overlapping bins $[E_{j},E_{j+1}]$ and defining $x_{j}=\int_{E_{j}}^{E_{j+1}}\bar{S}\left(E_{j}\right)dE$. Making the approximation that each $\mu_{s}$ is constant on each bin, we define the coefficient from the $j^{th}$ energy bin to the $i^{th}$ projection as $F_{i,j}=\exp\left\{-\sum_{s\in\Phi}\mu_{s}\left(E_{j}\right)L_{i,s}\right\}.$ Then (4) becomes

[TABLE]

Using this along with constraints to yield a normalized spectrum, the forward model is

[TABLE]

where $Y\in\mathbb{R}^{M}$ is a vector of $M$ normalized transmission measurements, $F\in\mathbb{R}^{M\times N_{e}}$ is a forward matrix over $N_{e}$ energy bins, $\tau$ is additive noise, and $x\in\mathbb{R}^{N_{e}}$ is the unknown vector discretization of $\bar{S}\left(E\right)$ that we seek to estimate.

2.2 Dictionary-based Model

Developing a dictionary-based method instead of directly estimating the spectrum at each energy bin has several motivations. When the discretization of $x$ is very fine, the projection model matrix $F$ has a large null space, leading to an underdetermined reconstruction, which requires a good initial estimate. Also, with a dictionary-based method, we can apply a sparsity-promoting penalty so that the dimensionality of the optimization problem can be significantly decreased and the time for reconstruction reduced.

More specifically, as indicated in (LABEL:equ:smp_model), we model $x$ as a discrete probability distribution, so that the entries of $x$ are nonnegative and sum to 1. We call this the simplex constraint and use $\mathcal{S}$ to represent the set of all such possible vectors.

We use a fixed, over-complete dictionary $D$ to represent the unknown normalized spectrum, $x$, as

[TABLE]

where $D$ is an $\mathbb{R}^{N_{e}\times N_{k}}$ matrix, each column $D_{*,k}\in\mathcal{S}$ represents a normalized basic spectrum. With this $D$, the transmission model can be rewritten as

[TABLE]

2.3 MAP Estimate

We use a Bayesian framework to estimate $\omega$ from transmission measurements under sparsity and simplex constraints. We define transmission weights using a diagonal matrix $\Lambda$, where we take $\Lambda_{i,i}=1/(y_{i}M)$. We then define the loss function $l(\omega)=\frac{1}{2}\|Y-FD\omega\|_{\Lambda}^{2}$, in which case the MAP estimate is given by

[TABLE]

where $\omega\in\mathcal{S}$ enforces $x\in\mathcal{S}$ since each column of $D$ is a simplex; the $L_{0}$ constraint ensures $\omega$ has no more than $C$ non-zero components.

2.4 Support Selection and Pairwise ICD update

To minimize the MAP cost function, we alternate between greedy support selection and a pairwise ICD update. Inspired by the Orthogonal Matching Pursuit (OMP)[16], our DictSE builds support set $\Omega=\{k:\omega_{k}\neq 0\}$ by adding one basic spectrum from the dictionary at a time and then updating the coefficients. However, the conventional support selection method in OMP does not account for simplex constraints on $\omega$. Further, the OMP method assumes that the dictionary atoms are normalized, whereas, in our spectral estimation problem, the product of forward matrix $F$ and the spectrum dictionary $D$ is not. Thus choosing a basis spectrum with conventional matching pursuit is inappropriate.

To describe our alternative method for support selection, we first define a function $l^{\prime}(k,\beta|\widehat{\omega})$ that measures the fit to data obtained by scaling the existing coefficients $\hat{\omega}$ by $\beta$ and using the remaining weight on the $k^{th}$ coefficient. That is,

[TABLE]

where $\beta\in[0,1)$ enforces $\omega\in\mathcal{S}$ and $\epsilon_{k}$ is a one-hot vector for the $k^{th}$ spectrum.

Then we select a new element from the dictionary by minimizing $l^{\prime}(k,\beta|\widehat{\omega})$ to obtain

[TABLE]

This minimization can be solved easily for each $k\notin\Omega$ since equation (10) is quadratic in the scalar $\beta$. In fact, defining $e_{k}=FD\widehat{\omega}-FD_{*,k}$, we have

[TABLE]

Using $\beta_{k}$ in equation (10) allows us to find $k^{*}$ that minimizes equation 10; this $k^{*}$ is then included in $\Omega$ with $\omega_{k^{*}}=1-\beta^{*}$ and the remaining $\omega_{k}$ scaled by $\beta^{*}$.

To rebalance the weights while enforcing the simplex constraint, we use pairwise ICD between the most recently added dictionary element with index $g$ and the remaining elements in $\Omega$, as shown in lines 6-12 of Algorithm 1. More precisely, after choosing $g$, we loop repeatedly over $k\in\Omega\setminus\{g\}$, in each case finding an optimal pairwise update $\omega\leftarrow\omega+\alpha^{*}(\epsilon_{k}-\epsilon_{g})$, where $\alpha^{*}$ is chosen by

[TABLE]

The constraints on $\alpha$ in the minimization ensure each $\omega_{k}\geq 0$. Since $l(\omega)$ is quadratic, $\alpha^{*}$ can be computed as below

[TABLE]

The pairwise ICD update will stop when the total update is less than $10^{-6}$. The algorithm is summarized in Algorithm 1.

3 Implementation

In this section, we describe the normalized transmission data $Y$, the projection matrix $F$, and the dictionary of spectra $D$, which are required to estimate the spectral response of an X-ray CT system using Algorithm 1.

3.1 Transmission data $Y$

We collected four CT datasets of different metal rods, $\Phi=\{Ti,V,Al,Mg\}$, at beamline 8.3.2 of the ALS. As shown in Fig. 2, for each dataset, we scanned a single rod. The same CT scanner setup was used to collect all datasets. Table 1 gives more detailed information about the X-ray CT measurements.

We scanned $15$ bright scans and $10$ dark scans and averaged them to obtain $\bar{I}_{\text{bright }}$ and $\bar{I}_{\text{dark }}$. Then, for each view, we normalized measurement data $I_{\text{scan }}$ to obtain $Y=\frac{I_{\text{scan }}-\bar{I}_{\text{dark }}}{\bar{I}_{\text{bright }}-\bar{I}_{\text{dark }}}$.

3.2 Forward Matrix $F$

For each dataset, the sample is a single rod that is solid and pure. To calculate $F_{i,j}$, we need the path length of the $i^{th}$ projection through this rod. To obtain the path length, we used filtered back projection to reconstruct the volume of the rod and then generated a mask to represent the object area. Using this mask, we calculated the path length for each projection and used this to calculate the forward matrix $F$.

3.3 Dictionary Generation $D$

As mentioned in equation (1), the spectral response can be modeled as the product of the X-ray source spectrum, filter response, and detector response, so we can create a dictionary by varying one or more of these elements. The X-ray source spectrum is fixed in this experiment; we used an estimate provided by the beamline scientists at the ALS Beamline 8.3.2. Therefore, in this implementation, we generated the dictionary by varying the filter and detector responses. Based on the spectrum models in Ref. [4], the filter response and detector response are determined by their material properties and thicknesses, which we vary as in Table 2 to generate our dictionary. By combining two groups of responses, we obtain an over-complete dictionary $D$ containing $60\times 36=2160$ normalized responses spectra.

4 Experimental Results

We compared our proposed DictSE method with a least-squares spectral estimation (LSSE) method provided by Livermore tomography tools (LTT) [17] on four datasets described in section 3.1. An initial spectrum for the LSSE method was generated by LTT using $3mm$ silicon as a source filter and $50\mu m$ $Lu_{3}Al_{5}O_{12}$ as a scintillator. We then evaluated the estimated spectra of DictSE and LSSE using leave-one-out cross-validation since we do not have a ground truth response.

Table 3 demonstrates that DictSE’s reconstructed spectra outperform the LSSE’s reconstructed spectra in NRMSE for all cross-validation cases. For each case $v\in\{1,2,3,4\}$, we computed NRMSE $=\frac{\|Y_{v}-\widehat{Y}_{v}\|_{2}}{\|Y_{v}\|_{2}}$ to compare transmission measurements $Y_{v}$ and transmission value of the forward model using estimated spectrum $\widehat{Y}_{v}=F_{v}D\widehat{\omega}_{v}$ on the validation rod.

Fig. 3 shows all four cases of cross-validation reconstructed spectra using both the DictSE and LSSE methods. For each case, DictSE’s estimated spectra are smoother than the LSSE’s. Also, from the shape of the reconstructed spectra over all cases, DictSE is less data-sensitive than LSSE.

5 Conclusion

This work provides a novel application of dictionary learning to X-Ray spectral estimation, allowing an efficient spectrum reconstruction from a vast dictionary obtained from CT datasets. Our method uses a greedy support selection method to do sparse coding followed by pairwise ICD to do minimization while enforcing a simplex constraint. Leave-one-out cross-validation experiments on four datasets demonstrated that our DictSE method outperforms the LSSE method in NRMSE.

6 ACKNOWLEDGMENTS

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 and LDRD project 22-ERD-011. The authors acknowledge Dula Parkinson for his support during beamtime. Charles Bouman was partially supported by the Showalter Trust, and Greg Buzzard was partially supported by NSF CCF-1763896.