1. Introduction

X-ray imaging is widely used for inspection in medical imaging, industrial testing, and homeland security. Conventional X-ray imaging systems, however, do not exploit the polychromatic spectrum of the X-ray source. Thus, the reconstructed gray-scale images are often insufficient for accurate medical diagnosis or material identification. Spectral X-ray computed tomography (CT) aims at revealing material characterization of an object by using information from the entire X-ray spectrum, which reveals the energy-dependent properties of the linear attenuation coefficients [1]. Various methods can be used to obtain information from different X-ray energies. Dual-energy CT systems obtain two data-sets with distinct attenuation characteristics by using two sequential CT scans with different X-ray tube voltages [2], or a set of sandwiched detectors capable of discriminating between low and high X-ray energies [3]. Nonetheless, these systems are inaccurate at low energies or for scenes containing various materials with characteristic K-edges, as more than two spectral measurements are needed [4, 5]. Photon counting detectors are able to obtain attenuation data from more than two X-ray energy bins, therefore allowing more materials to be distinguished in a single acquisition [6]. Yet, such detectors are prohibitively costly for some applications, and the measured data have low signal to noise ratio (SNR) due to the detector’s narrow bin bandwidth and quantum noise which has limited their application in a clinical setting [7–9]. Furthermore, the complexity of the sensor layers limits the resolution of the detectors, which is compromised when more energy bins are needed [10]. Integrating detectors typically offer higher pixel count but suffer from higher dark noise and lower quantum efficiency (QE). Newer CMOS integrating detectors, however, are making significant progress in resolution and high-emission efficiency [11] to improve the overall QE. Rakvongthai et al. proposed a cost-effective spectral CT system consisting of K-edge filters integrated with a conventional X-ray system [12] based on the Ross filter method [13]. A K-edge filter is a material consisting of a high-Z element, i.e. an element with a high atomic number, such as tantalum, tungsten, or molybdenum, which sharply cuts off part of the X-ray spectrum above the element’s K-shell electrons’ binding energy. Thus, using different K-edge filters in sequential scans, mono-energetic X-ray flux data can be obtained by a combination of the sequential filtered measurements (see Fig. 1). The balanced filter method proposed by Rakvongthai et al. [12] uses this principle to reconstruct energy binned images overcoming the limitations of photon counting detectors, at the expense of longer acquisition times and higher radiation dose given that a complete scan is needed for each K-edge filter used. This sequential acquisition method has been widely used in X-ray crystallography [14], spectroscopy [15] and more recently in X-ray CT monochromatization in medical imaging [16].

 

Fig. 1 (a) 80 keV X-ray source spectrum (blue), filtered with Ce (red) and Dy (cyan). Dy and Ce are a Ross filter pair; thus, a quasi-monochromatic spectrum can be obtained by subtracting the filtered spectrum corresponding to the lower K-edge filter, Ce, from the spectrum, obtained using the higher K-edge filter, Dy as shown in (b). Vertical lines in the spectra indicate the K-edges of Ce and Dy, 40.4 and 53.8 keV, respectively.

Download Full Size | PDF

 

Fig. 2 Top: Subsampling strategies for CSXI. (a) Uniform view (UV-CSXI): K-edge filters are assigned sequentially every F view angles (F is the number of K-edge filters used); (b) Random view (RV-CSXI): A K-edge filter is assigned randomly to every view angle; (c) Uniform detector (UD-CSXI): the same pattern of K-edge filters repeats every F coded aperture elements; (d) Random detector (RD-CSXI): random coded apertures with multiple K-edge filters are used at each view angle. Bottom: Spectrally coded sinograms, color-coded according to the material used to filter the X-ray that impinges on a particular detector element. Zoomed area is shown for each sinogram, depicting the multiplexing structure.

Download Full Size | PDF

The work in this paper presents a different approach for spectral CT which relies on spatially and spectrally coded illumination of objects, yet measures the projections using standard X-ray integrating detectors during a single scan. Since the multiple energy bins are coded and compressively measured at the integrating detector, the new approach is coined compressive spectral X-ray imaging (CSXI). Various spectral and spatial coding strategies for structured illumination are explored. Subsamplingin CSXI is accomplished by either view angle spectral subsampling, spatial subsampling enabled by coded apertures placed at the source or detector side, or both. In all of these strategies, the X-ray source spectrum is spectrally coded using balanced Ross filters while the spatial information is coded by placing the K-edge filters in a way that the view angles or detectors are subsampled. Reducing the number of view angles is equivalent to placing a different K-edge filter at each X-ray source position; the filter placement can be uniform or random. On the other hand, subsampling of the detectors can be achieved by using standard block-unblock or K-edge coded apertures placed at the source or detector side. Figures 2 and 3 depict the various coding strategies in CSXI. Figure 2(a) shows spectral view subsampling by sequentially switching a set of spectral filters along each view. Figure 3(a) depicts the same strategy but in addition, it places a block-unblock coded aperture to subsample the detectors using structured illumination. Figures 2(b) and 3(b) depict the same strategies as those in Figs. 2(a) and 3(a), except that the spectral filter placement is random rather than sequential. Figure 2(c) shows a strategy where spectral and spatial subsampling is attained by placing a uniformly distributed K-edge coded aperture in front of the source. Figure 3(c) shows the same approach with the addition of blocking elements on the coded aperture for additional detector subsampling. Finally, Figs. 2(d) and 3(d) illustrate the same strategy except that the elements of the K-edge coded apertures and the blocking elements, are all placed randomly. The eight types of CSXI configurations depicted in Figures 2 and 3 are referred to as Uniform View (UV-CSXI), Random View (RV-CSXI), Uniform Detector (UD-CSXI), Random Detector (RD-CSXI), Uniform View Block/Unblock (UVB-CSXI), Random View Block/Unblock (RVB-CSXI), Uniform Detector Block/Unblock (UDB-CSXI), and Random Detector Block/Unblock (RDB-CSXI), respectively.

Variations of the above approaches are possible which can be useful in simplifying the implementation. For instance, a single coded aperture can be translated as the angle is rotated. Another approach can use a set of coded apertures that are placed on the other side of the object, just on top of the detector. Most importantly, in all of the above CSXI sampling strategies, a single rotation scan is used and the measurements are captured with a conventional integrating detector, circumventing the limitations of other spectral tomography approaches — albeit the simpler scanning mechanism is compensated by added computational complexity. The advantages of the proposed CSXI approaches lie in: (i) structured illumination of the object leading to reduced radiation exposure for patients in medical scans; (ii) spectral and spatial subsampling which can overcome or ameliorate limited angle view constraints in security inspection applications; (iii) subsampling which does not lead to image quality degradation since the principles of compressive sensing are used based on the understanding that real objects have sparse representations in some set of basis functions, which indicates that objects can be defined in a smaller subspace than their number of pixels or voxels.

 

Fig. 3 Top: X-ray downsampling in CSXI with a subsampling rate of 2. A random block/unblock coded aperture is assigned to each view angle in the (a) UVB-CSXI and (b) RVB-CSXI configurations. And, blocking elements are introduced in the K-edge coded apertures to subsample the measurements for the (c) UDB-CSXI and (d) RDB-CSXI configurations. Bottom: The spectrally coded sinograms in each case are color-coded according to the material used to filter the X-ray that impinges on a particular detector element. A zoomed area of each sinogram is shown to illustrate their multiplexing structure.

Download Full Size | PDF

The effective linear attenuation coefficients can be reconstructed from the set of measurements corresponding to a particular filter by using compressive sensing algorithms to solve the monochromatic inverse problem as shown in [17]. Yet, these reconstructions do not provide the spectral information of the object. Furthermore, the forward model in CSXI cannot be written as a linear system of equations in terms of the energy binned linear attenuation coefficients since the measurements from the integrating detectors correspond to a nonlinear mapping of the spectral datacube. Thus, conventional compressed sensing reconstruction algorithms cannot be directly applied to reconstruct the energy binned linear attenuation coefficients from the filtered measurements. Note, however, that each CSXI ray is associated with a particular K-edge filter. Having a set of F different filters in the coded aperture, the filtered measurements obtained at the detector provide F subsampled multiplexed sinograms. Hence, the goal is to develop efficient computational methods to recover the points that are not observed. Then, the physics governing the spectral X-ray projection measurements together with compressive sensing are used to estimate the energy binned sinograms. Subsequently, the alternating direction method of multipliers (ADMM) is used to solve the inverse problem by exploiting the inherent sparsity of X-ray images in the spatial domain and the low-rank structure of the spectral X-ray datacube [18]. The processed image can be further decomposed into basis functions to obtain material based images, which are useful for diagnosis and analysis in medical and/or security applications. We compare the proposed sampling strategies and evaluate their performance at different subsampling rates. Simulations show that CSXI reconstructions exhibit up to 7 dB gain when compared to the reconstructions attained by the sequential balanced filter framework in [12] when using an equivalent radiation dose. The generalization to 4D spectral tomography based on cone-beam projections, the possible resolution miss-match between the detector and the coded apertures, and the optimization of the various coded apertures used in CSXI are topics of current research. The remainder of this paper is organized as follows: Section 2 describes the general sampling problem and presents the various sampling strategies for CSXI. The CSXI forward model is developed in Section 3, as well as, the formulation of the inverse problem and the algorithms to reconstruct the spectrally resolved images. Section 4 presents simulation results for a fan-beam scenario for all the subsampling strategies and Section 5 presents the conclusions and future work.

2. Problem statement

2.1. Preliminaries

When a polychromatic X-ray beam passes through an object, the X-ray attenuation at each energy, E, is obtained as $I\left(E\right)=S\left(E\right) Q\left(E\right) P\left(E\right) \exp (-\underset{\mathcal{l}}{{\displaystyle \int }}\mu \left(\mathcal{l},E\right)d\mathcal{l})$ [19], where $\mu \left(\mathcal{l},E\right)$ is the linear attenuation coefficient of the object at position $\mathcal{l}$, $S\left(E\right)$ is the X-ray source spectrum, $Q\left(E\right)$ is the energy-dependent detector response, and $P\left(E\right)$ is the effective energy flux. For photon counting detectors $P\left(E\right)=1$ and for conventional integrating detectors $P\left(E\right)=E$. That is, photon counting detectors reflect the number of photons in each energy bin; while integrating detectors accumulate the intensity over the entire X-ray energy spectrum and thus the contribution of each photon to the reading is weighted by its original energy [10]. Let $I_0\left(E\right)=S\left(E\right)\cdot Q\left(E\right)\cdot P\left(E\right)$ be the intensity measured at energy E without any object in front of the X-ray source; then, for a conventional X-ray imaging system, with integrating detectors, the intensity registered at the $j^{th}$ detector element is given by:

If P is the total number of X-ray source positions and M is the number of detector elements in the detector array, then $j=1,\dots ,MP$. The integrated gray-scale measurements in Eq. (1) cannot be used to reconstruct directly the energy binned images $\mu \left(\mathcal{l},E\right)$, as they do not provide spectrally resolved data [20]. Photon counting detectors and monochromatic synchrotron X-ray sources are able to provide such energy selective measurements; however, these devices are costly. As mentioned previously, a cost-effective alternative to obtain quasi-monochromatic X-ray spectra are Ross filters [16, 21]. In principle, Ross filters produce a monochromatic effect by using two balanced K-edge filters that consist of materials with nearly adjacent atomic numbers whose thicknesses are carefully matched such that the transmitted spectra are identical for all photon energies except in the narrow energy bin between their respective K-edges [22, 23]. The quasi-monochromatic measurements are then obtained by subtracting the X-ray intensity acquired using the filter with the lower K-edge from that with the higher K-edge in the Ross pair [16]. The bandwidth of the energy bin is thus defined by the difference between the K-edge energies of the two Ross pair elements. Figure 1(a) illustrates the process to obtain a quasi-monochromatic X-ray spectrum using the Ross pair constituted by Dysprosium (Dy) and Cerium (Ce); by subtracting the spectrum filtered by Ce from the spectrum filtered by Dy a quasi-monochromatic spectrum between 40.4 keV and 53.8 keV is obtained (Fig. 1(b)). Rakvongthai et al. proposed a multiple Ross filter approach for spectral CT which uses a conventional X-ray fan-beam CT system equipped with five K-edge filters, that are used separately to obtain five sequential measurement sets [12]. To compensate for the unbalanced material thicknesses, the authors use the inverse of a transmission matrix that models the spectral response of the filters in each energy bin. Finally, the linear attenuation coefficients for each energy bin are obtained by applying the filtered back-projection (FBP) algorithm to the respective energy binned sinogram. This approach provides a cost-effective alternative to photon counting detectors, but the acquisition time, and thus the radiation dose increase since one complete scan per filter is required. In general, to obtain measurements corresponding to K different energies, a set of $K+1$ sequential filtered projections are needed.

2.2. Spectral CT subsampling strategies

Subsampling is desired in spectral CT to reduce both the scanning time and the radiation dose while maintaining the reconstruction quality. In monochromatic X-ray imaging, Kaganovsky et al. introduced subsampling strategies for monochromatic X-rays in medical CT scanner geometries [24]. These systems require coded apertures to spatially block the X-ray projections in some fashion or subsample the number of source positions to obtain projections at a lower number of angles. In contrast, this work explores subsampling strategies not only in the spatial dimension as done in [24] but also in the X-ray source spectrum. K-edge filters, block-unblock coded apertures and/or K-edge coded apertures are placed in a way that the source positions or detectors are subsampled and the X-ray source spectrum is modified. Figure 2 depicts the sketches of the four subsampling configurations proposed for CSXI as well as their corresponding intensity sinograms for F = 5 filters. The four subsampling strategies defined previously and depicted in Figs. 2 and 3 are: (a) Uniform view (UV-CSXI), in which a particular K-edge filter is assigned every $F^{th}$ angle, such that the set of K-edge filters are uniformly distributed in all X-ray source positions; (b) Random view (RV-CSXI), in which a particular filter is assigned randomly to every view angle; (c) Uniform detector (UD-CSXI), in which the same pattern of K-edge filters repeats every F coded aperture elements; and (d) Random detector (RD-CSXI), in which random coded apertures with multiple random K-edge filters are used to obtain spatially and spectrally coded illumination projections of the object. At the bottom of Fig. 2, the sinograms for each subsampling strategy are shown. Each column of the sinogram represents a view angle of the system and each row represents the detector element location in the detector array. A different color is used for each filter for illustration purposes. Note, that each configuration has a very characteristic sampling pattern. The sinograms associated with the strategies that do not use coded apertures, i.e. UV-CSXI and RV-CSXI, have a particular color assigned to each column which corresponds to each of the view angles Fig. 2(a)-2(b). On the other hand, in the sinograms associated with the strategies that use K-edge coded apertures, i.e. UD-CSXI and RD-CSXI, each pixel corresponds to a different filter as depicted in Fig. 2(c)-2(d).

Further X-ray subsampling can be achieved in two ways: reducing the number of view angles or reducing the number of measurements at the detector array by using block-unblock coded apertures in front of the X-ray source to create structured X-ray bundles. As will be shown in Section 4, random detector subsampling in each view angle presents the reconstructions with the lowest error. To this end, block/unblock coded aperture structures are placed after the K-edge filter in the UVB-CSXI and RVB-CSXI architectures as shown in Figs. 3(a)-3(b), respectively. For the UDB-CSXI and RDB-CSXI configurations, the reduction of measurements is achieved by adding blocking elements to the K-edge coded apertures as shown in Figs. 3(c)-3(d), respectively. The corresponding sinograms are shown in the bottom row of Fig. 3.

3. CSXI forward model and the inverse problem

All of the CSXI configurations can be seen as a particular case of RD-CSXI; specifically, UV-CSXI corresponds to having the same K-edge filter in all the K-edge coded aperture elements every $F^{th}$ view angle; RV-CSXI corresponds to assigning the same K-edge filter to all K-edge coded aperture elements of a particular view angle; and UD-CSXI corresponds to repeating the K-edge coded aperture pattern every F coded aperture elements. Thus, without loss of generality, we consider next the general CSXI forward imaging model where K-edge coded aperture patterns are placed in front of the X-ray source at each position s_p, for $p=1,\dots ,P$, as shown in Figs. 2(d) and 3(d). When a filter f is placed in front of an X-ray source, the intensity of an X-ray beam at the energy E after passing through the filter is given by $I_f\left(E\right)=I_0\left(E\right)\exp \left[-{\mu }_f\left(E\right){\delta }_{j,f}\right]$, where ${\mu }_f\left(E\right)$ is the linear attenuation coefficient of the filter f at the energy E, $I_0\left(E\right)$ is the X-ray source intensity at energy E, and ${\delta }_{j,f}$ is the length of the intersection of the $j^{th}$ X-ray beam with the filter f. The latter is given by ${\delta }_{j,f}={\rho }_f/\cos ({\psi }_j)$, where ρ_f is the thickness of the filter and ψ_j is the angle between the normal of the filter and the $j^{th}$ X-ray beam. Thus, the filtered measurement at the $j^{th}$ integrating detector element is given by:

To formulate the inverse problem, the linear attenuation coefficients are discretized and then, as shown in Appendix A, Eq.(2) can be rewritten as

This minimization problem is non-linear and the set of measurements $I_j^f$ are highly undersampled, therefore high quality reconstructions cannot be obtained using a conventional gradient minimization approach. This paper presents a method to estimate the energy binned sinograms from the integrating detector measurements to then reconstruct the energy binned linear attenuation coefficients ${\mathbf{u}}^k$. Note that the set of filtered projections measured by the integrating detector, in essence, builds up a set of F tiled subsampled sinograms ${\mathbf{I}}^f$, as shown in Fig. 4(a). The inverse imaging problem would aim at sinogram inpainting to recover all the missing points in the log-normalized version of ${\mathbf{I}}^f$, defined as ${\mathbf{y}}^f$, for $f=1,\dots ,F$ (Fig. 4(b)). Unfortunately, ${\mathbf{I}}^f$ and ${\mathbf{y}}^f$ are both highly undersampled such that straight inpainting leads to poor image reconstruction. A better approach is needed, one that exploits the physics behind the spectral sensing together with the low-rank and sparse representation of the data. This is accomplished as follows. First, the sinogram inpainting problem is solved for each incomplete log-transformed sinogram ${\mathbf{y}}^f$ which yields the effective linear attenuation coefficients using filter f, ${\mathbf{x}}^f$ (Fig. 4(c)), and after the inverse log-normalization operation the estimate ${\widehat{\mathbf{I}}}^f$ is obtained (Fig. 4(d)). These procedures are detailed in Section 3.1. The next phase, once the filtered intensities ${\widehat{\mathbf{I}}}^f$ have been estimated at all points, is to take this set of measurements into the energy domain using the transmission matrix B, as shown in Fig. 5, to obtain the estimated log-transformed energy binned sinograms ${\mathbf{z}}^k$. These steps will be detailed in Section 3.2. Finally, the effective linear attenuation coefficients at each energy bin ${\mathbf{u}}^k$ can be reconstructed from these sinograms using a joint sparse and low-rank optimization method as detailed in Section 3.3. To show the inverse problem clearly, each phase of the reconstruction steps is described in separate subsections.

 

Fig. 4 Filtered sinograms inpainting. (a) I^f, measured undersampled sinograms; (b) y^f , log-normalized undersampled sinograms; (c) x^f , reconstructed effective linear attenuation coefficients; (d) ${\widehat{\mathbf{I}}}^f$, final estimated intensity sinograms, for $f=1,\dots ,F$.

Download Full Size | PDF

3.1. Filtered sinograms inpainting

Each measurement $I_j^f$ is associated with a particular filter f; for example, Fig. 4(a) illustrates the spectrally coded sinograms when F = 5 filters are used. If F different filters are used in the coded apertures, the CSXI measurements lead to F distinct under sampled sinograms ${\mathbf{I}}^f$, which are then normalized to obtain the log-normalized sinograms ${\mathbf{y}}^f$ as shown in Fig. 4(b). Then, each subsampled set of measurements ${\mathbf{y}}^f$ is used to reconstruct the effective linear attenuation coefficients ${\mathbf{x}}^f$ using regularization constraints to account for the missing data as depicted in Fig. 4(c). Subsequently, the effective linear attenuation coefficients ${\mathbf{x}}^f$ can be re-projected to the sinogram space using the information of the CT system matrix H to obtain the missing points in ${\mathbf{I}}^f$ as illustrated in Fig. 4(d). This process is described in more detail as follows. Recall that an effective X-ray energy and an effective linear attenuation coefficient can be defined to treat the polychromatic X-ray problem using the monochromatic X-ray model [19]: $I=I_0\exp (-\underset{\mathcal{l}}{{\displaystyle \int }}\overline{\mu }\left(\mathcal{l}\right)d\mathcal{l})$, where I₀ is the intensity of the X-ray source at the effective energy, I is the measured intensity on the detector, and $\overline{\mu }\left(\mathcal{l}\right)$ are the effective linear attenuation coefficient along the projected X-ray through the object. Hence, Eq. (2) can be re-written for each filter f as:

While many algorithms exist to solve Eq. (6), the Two-Step Iterative Shrinkage/Thresholding (TwIST) algorithm [26] is used here to reconstruct ${\mathbf{x}}^f$ and the Total Variation (TV) norm of ${\mathbf{x}}^f$ is used as the regularization term. Then, the inpainted sinogram ${\widehat{\mathbf{y}}}^f$ is attained using the complete system matrix to estimate the missing points, i.e. ${\widehat{\mathbf{y}}}^f=\mathbf{H}{\mathbf{x}}^f$. The estimated intensity at the detector element j using the $f^{th}$ filter is then obtained as ${\widehat{I}}_j^f={\mathcal{N}}^{-1}\left({\widehat{y}}_j^f\right)=I_0^f\exp (-{\widehat{y}}_j^f)$. Note, this process might change the values of ${\mathbf{I}}^f$ at the locations that had been originally measured. Therefore, the original measurements are substituted into the estimated intensity sinograms; mathematically, ${\widehat{\mathbf{I}}}_{{\mathrm{\Omega }}_f}^f={\mathbf{I}}^f$, where Ω_f corresponds to the set of values of j where the filter f was used.

 

Fig. 5 Monochromatic sinogram estimation process. (a) ${\widehat{\mathbf{I}}}^{\mathrm{f}}$, stacked estimated intensity measurement with column vectors r_j ; (b) J, energy binned intensity estimates with column vectors d_j; (c) z^k, log-transformed energy binned sinograms.

Download Full Size | PDF

3.2. Energy binned sinograms estimation

Given the estimated sinograms ${\widehat{\mathbf{I}}}^f$ (Fig. 5(a)), in principle, one could find the monochromatic intensities using the balanced Ross filter method at each detector element j by subtracting the intensities corresponding to a Ross filter pair. As mentioned before, however, given the error introduced by the estimation process described in Section 3.1, subtracting the measurements to obtain monochromatic sinograms leads to poor reconstructions. Instead, the physical phenomenon in the sensing of X-rays is used to attain a more accurate estimation of the energy binned measurements.

Given a set of Ross filters, based on their K-edges, the X-ray spectrum is divided into $F+1$ energy bins $\left(0,E_1\right],\left(E_1,E_2\right],\dots ,\left(E_K,E_F\right],\left(E_F,E_T\right]$, where E_f for $f=1,\dots ,F$ is the K-edge energy of the filter f, $K=F-1$, and E_T is the maximum energy of the X-ray source. Figure 6 depicts an example for F = 3 filters, Cerium (Ce), Dysprosium (Dy), and Erbium (Er). Note the spectrum is sharply filtered after the K-edge energy of each element, 40.4 keV (E₁) for Ce, 53.8 keV (E₂) for Dy, and 57.5 keV (E_F) for Er. Furthermore, recall the thicknesses of the filter pairs are calculated such that the transmitted spectra is almost identical except in the energy band between their K-edges. Thus, for the filter pair Dy-Ce the filtered spectra between $\left(0,E_1\right]$ and $\left(E_2,E_T\right]$ is almost identical as shown in Figs. 6(b) and 6(c); and for the filter pair Er-Dy, the filtered spectra between $\left(0,E_2\right]$ and $\left(E_F,E_T\right]$ is almost identical as shown in Figs. 6(c) and 6(d). Thus, for the set of filters in the example, the spectra obtained with all the filters in the bands $\left(0,E_1\right]$ and $\left(E_F,E_T\right]$ is almost identical. Generalizing for a set of F K-edge filters with balanced thicknesses, the energy spectra in the intervals $\left(0,E_1\right]$ and $\left(E_F,E_T\right]$ are almost identical for all the filters.

The filtered intensities in the integrating detector are a linear combination of the desired monochromatic intensities and can be written as $I_j^f={\displaystyle \sum }_{k=1}^{F+1}{\mathbf{B}}_j\left(f,k\right)d_j^k$, where ${\mathbf{B}}_j\left(f,k\right)=\left(\underset{E_k}{\overset{E_{k+1}}{{\displaystyle \int }}}\exp \left(-{\mu }_f\left(E\right){\delta }_{j,f}\right)dE\right){\left(E_{k+1}-E_k\right)}^{-1}$, and $d_j^k$ is the monochromatic intensity at the detector element j at energy bin k, $\left(E_k,E_{k+1}\right]$ (see Appendix A). Given that the energy spectrum in the intervals $\left(0,E_1\right]$ and $\left(E_F,E_T\right]$ are very similar in value for all the filters, the recovery of these bands (k = 0 and $k=F$) is not possible. Let ${\mathbf{r}}_j={\left[{\widehat{I}}_j^1,{\widehat{I}}_j^2,\dots ,{\widehat{I}}_j^F\right]}^T$ be the estimated filtered intensities, and ${\mathbf{d}}_j={\left[d_j^0,d_j^1,d_j^2,\dots ,d_j^K,d_j^F\right]}^T$ be the energy binned intensities at the $j^{th}$ detector element. Then, the filtered intensities at the $j^{th}$ detector can be expressed as ${\mathbf{r}}_j={\mathbf{B}}_j{\mathbf{d}}_j$, where ${\mathbf{B}}_j\in {ℝ}^{F\times \left(F+1\right)}$ is a weighting matrix that accounts for the spectral response of the filters at each energy bin. It should be noted that only one filtered intensity is measured by the CSXI system at the $j^{th}$ detector, the remaining $F-1$ filtered intensities are estimates obtained by the process described in Section 3.1 and may contain errors. Thus, the pseudo-inverse of the matrix ${\mathbf{B}}_j$ cannot be used directly to obtain the monochromatic intensities. In order to account for the noisy estimates, ${\mathbf{r}}_j$, a joint estimation of the monochromatic X-ray intensity sinograms using a gradient descent algorithm is used to solve the inverse problem:

 

Fig. 6 (a) 80 keV X-ray source spectrum, filtered with (b) Ce, (c) Dy, and (d) Er. Ce and Dy, and Dy and Er are Ross filter pairs. Note the filtered spectra for all filters has the same value in the bins $\left[0,E_1\right)$ and $\left(E_f,E_T\right)$ highlighted in black for all the figures.

Download Full Size | PDF

3.3. Reconstruction of the energy binned images

The monochromatic intensity sinograms depicted in Fig. 5(b), can be used to reconstruct the energy binned linear attenuation coefficients $\mu \left(\mathcal{l},k\right)$ at energy bin k, for $k=1,\dots ,K$. Using the monochromatic X-ray Beer-Lambert law, the energy binned intensities at the $j^{th}$ detector can be written as ${\widehat{d}}_j^k=d_{j0}^k \exp \left(-\underset{\mathcal{l}}{{\displaystyle \int }}\mu \left(\mathcal{l},k\right)d\mathcal{l}\right)$, where $d_{j0}^k$ is the quasi-monochromatic intensity at energy k. Normalizing and discretizing the line integrals and the linear attenuation coefficients $\mu \left(\mathcal{l},k\right)$, the energy binned measurements ${\mathbf{z}}^k={[z_1^k,z_2^k,\dots ,z_{MP}^k]}^T$ (Fig. 5(c)) are given by the following linear model: ${\mathbf{z}}^k=\mathbf{H}{\mathbf{u}}^k$, where ${\mathbf{u}}^k\in {ℝ}^{N^2\times 1}$ is the vectorized linear attenuation coefficient $\mu \left(\mathcal{l},k\right)$ of the object at the $k^{th}$ energy bin and H is the CT system matrix accounting for the hardware settings. Recall the columns of the system matrix are associated with the pixels of the object while its rows correspond to the detector elements. For the full sampling scenarios (Fig. 2) H contains the information of all the detectors. For the subsampling scenarios depicted in Fig. 3, the rows of H correspond to the detectors that are not paired with a blocking element.

Given the error introduced by the multiple sinogram estimation phases depicted in Figs. 4-5 and the spatial subsampling performed in the different strategies (Fig. 3), regularization constraints based on the structure of the data are used to reconstruct ${\mathbf{u}}^k$. Here, a joint sparse and low-rank optimization method is used, but other approaches such as [27], which seeks to minimize the nuclear norm in patches of the spectral images can be used as well. The approach used here jointly minimizes the ${\mathcal{l}}_1$ norm of the sparse representation of the datacube, and its nuclear norm. Furthermore, the tensor modeling previously proposed by Gao et al. in [10] and [28] is adapted to formulate the reconstruction problem as follows:

Vectorization: Let $\text{vec}\left(\mathcal{U}\right)$, be the operator $\text{vec}:{ℝ}^{Q_1\times Q_2\times \dots \times Q_N}\to {ℝ}^{Q_1\cdot Q_2\dots Q_N}$ that stacks the entries of a tensor in reverse lexicographical order into a $Q_1\cdot Q_2\dots Q_N$-long column vector [29].

Unfolding: For a tensor $\mathcal{U}$, the matrix unfolding denoted by ${\mathbf{U}}_{\left(n\right)}={\text{Unfold}}_{n-th}\left(\mathcal{U}\right)$ is defined by the mapping of: ${ℝ}^{Q_1\times Q_2\times \dots \times Q_N}\to {ℝ}^{Q_n\times \left(Q_{n+1}\cdot Q_{n+2}\dots Q_N\cdot Q_1\cdot Q_2\dots Q_{n-1}\right)}$. The resulting matrix ${\mathbf{U}}_{\left(n\right)}$ is a 2D matrix with Q_n rows and $Q_{n+1}\cdot Q_{n+2}\dots Q_N\cdot Q_1\cdot Q_2\dots Q_{n-1}$ columns [30].

Folding: The operator $\mathcal{U}={\text{Fold}}_{n-th}\left({\mathbf{X}}_{\left(n\right)}\right)$ is the inverse operator of unfolding [10].

Using these concepts, the fan beam model used here can be generalized to multiple geometries by representing the spectral linear attenuation coefficients with a tensor $\mathcal{U}\in {ℝ}^{Q_1\times Q_2\times \dots \times Q_{\xi +1}}$, where ξ is the dimensionality of a single energy X-ray image. For instance, ξ = 2 for fan-beam and parallel beam CT, ξ = 3 for cone-beam CT. Without loss of generality, this paper treats the inverse problem for a fan-beam system; therefore, the tensor $\mathcal{U}\in {ℝ}^{N\times N\times K}$ is a multi-dimensional matrix, where the first two dimensions are spatial, the third is the energy dimension, and $\text{vec}\left(\mathcal{U}\right)={\left[{\mathbf{u}}^{1\mathrm{T}}|{\mathbf{u}}^{2\mathrm{T}}|\dots |{\mathbf{u}}^{K\mathrm{T}}\right]}^T$. The spectral CT model can thus be generalized as:

3.3.1. Energy binned image reconstruction via ADMM

An additional regularization constraint is added to the inverse problem to take into account the correlation across spectral channels. This is achieved by introducing a low-rank constraint on $\mathcal{U}$ [10, 18] as:

The problem is then split into three different sub-problems summarized in Algorithm 1. Step 2 of the algorithm is an ${\mathcal{l}}_2$ minimization problem which is solved using the conjugate gradient (CG) algorithm with a fixed step size, and as introduced in [18] the solution for Step 3 and Step 4 are respectively given by:

Algorithm 1. Image reconstruction via the ADMM

View Table | View all tables in this article

3.3.2. Initialization

The X-ray attenuation coefficients vary smoothly and have sharp boundaries. TV regularization is known to be a good reconstruction method for such cases, as it is robust to noise and the edges in the reconstructions are well-defined [32]. Hence, in order to find a suitable initialization point for the ADMM algorithm, the energy binned images ${\mathbf{u}}^k$ are reconstructed from the corresponding quasi-monochromatic sinograms by solving the following minimization problem independently for each energy bin:

 

Fig. 7 (a) Forbild thorax phantom [33]. (b) Energy spectra of an 80 keV X-ray source seen through the five different filters used in the CSXI simulations [Molybdenum (Mo), Cerium (Ce), Dysprosium (Dy), Erbium (Er) and Tungsten (W)]. (c) Quasi-monochromatic energy binned spectra obtained subtracting the filtered spectra of the corresponding Ross filter pairs Ce-Mo, Dy-Ce, Er-Dy, and W-Er.

Download Full Size | PDF

Table 1. Ross filter pairs and respective energy bins

View Table | View all tables in this article

4. Simulation results

ACT configuration with a flat 1D detector strip with $M=512$ elements, a fan beam X-ray source, and $P=512$ view angles was used in the simulations. The scanned object was a 256 × 256 modified Forbild thorax phantom [33] generated using the CONRAD software [34]. The modified phantom shown in Fig. 7(a) consists of eight different materials to simulate (1) air, (2) lung, (3) average soft tissue, (4) heart (blood), (5) iodine$+$ blood, (6) artery (blood), (7) bone, and (8) marrow. The mass attenuation coefficients were obtained from the National Institute of Standards and Technology (NIST) X-ray attenuation databases available in [35]. The distance from the center of rotation to the X-ray source as well as to the detector strip is set to 40 cm and the detector length is set to $41.3$ cm with $M=P=2N=512$, to match the full sampling scenario described by Jorgensen et al in [36]. Five different filter materials [Molybdenum (Mo), Cerium (Ce), Dysprosium (Dy), Erbium (Er) and Tungsten (W)] are used as K-edge filters, i.e. F = 5. After matching the Ross filter pairs to obtain energy binned sinograms as described in [16], the thickness of the Mo, Ce, Dy, Er, and W filters were set to be 74.7, 52.8, 30.6, 26.7 and 9.9 μm, respectively. The resulting energy bins were Ce-Mo [20.0-40.4], Dy-Ce [40.4-53.8], Er-Dy [53.8-57.5] and W-Er [57.5-69.5] keV, for $k=1,2,3,4$ respectively, as detailed in Table 1. The X-ray filtered spectra $I_f\left(E\right)$ were simulated at 80 keV using the Spektr software [37] and a copper filter of 0.25 mm was placed in front of the X-ray source to filter out the energies lower than 20 kV. Figure 7(b) depicts the energy spectra of the simulated 80 kV X-ray source seen through the five different filters and Fig. 7(c) depicts the spectra of the energy bins obtained using the Ross filter pairs.

 

Fig. 8 (a) Effective linear attenuation coefficient images at energy bin k = 3. Reconstructions using a subsampling rate of 8 for (b) UV-CSXI when view angles are subsampled and (c) UVB-CSXI (block-unblock coded apertures are used to subsample the detectors). The PSNR and zoomed versions of the highlighted regions are shown for each reconstruction. (d-e) Normalized absolute errors for (b) and (c), respectively.

Download Full Size | PDF

The projection data were simulated using the ASTRA tomography toolbox [38] for the setup previously described and energy weighted integrals over 1keV spectral steps were obtained for each filtered measurement according to Eq. (2). For simplicity, the coded apertures are assumed to have the same number of elements as the detector array and the coded aperture pitch is fixed to obtain a one to one correspondence with the detector elements, for both the K-edge coded apertures and the monochromatic coded apertures. However, when the coded aperture elements do not have a one to one correspondence with the detector array, computational methods can be implemented to address the pixel mismatch. A possible approach is to adapt the miss-match algorithm proposed in [39]. Another approach is to model the system in which the coded aperture mask is located after the object and close to the detector which can be shown to have the same forward model.

Table 2. CSXI PSNR at each k energy bin for the Thorax Forbild Phantom

View Table | View all tables in this article

View angle subsampling is often used to reduce the number of measurements in X-ray imaging systems due to its straightforward implementation. Simulations to compare the performance of the UV-CSXI configuration for uniform view subsampling (no coded apertures) and detector subsampling UVB-CSXI (block-unblock coded apertures as shown in Fig. 3(a)) are performed. In the ADMM algorithm the energy binned CT image is represented on the sparse basis $\mathbf{\Psi }={\mathbf{\Psi }}_{DCT}\otimes {\mathbf{\Psi }}_W$, where $\otimes$ is the Kronecker product, ${\mathbf{\Psi }}_{DCT}$ is the discrete cosine transform (DCT) basis in the energy domain, and Ψ_W is the 2D Haar wavelet basis in the spatial domain. The reconstructions are compared to the effective linear attenuation coefficients at each energy bin calculated using the NIST tables [35]. The peak signal to noise ratio (PSNR) is used to evaluate the reconstructions since it is suitable for comparing restoration results [17]. Figure 8(a) depicts the effective linear attenuation coefficients of the thorax phantom for k = 3 and Figs. 8(b) and 8(c) depict the reconstructions obtained for a subsampling rate of 8 for UV-CSXI when the number of view angles is subsampled and UVB-CSXI, respectively. Both systems measure the same number of line integrals, however, as depicted in Fig. 8, UVB-CSXI is significantly better than UV-CSXI when the view angles are subsampled. Particularly, in the zoomed versions and the normalized absolute error of the reconstructions (Figs. 8(c)-8(d)) it can be seen that the objects inside the phantom cannot be accurately reconstructed when the view angles are subsampled. It can be concluded that detector subsampling leads to higher reconstruction quality in CSXI. Hence, for the remaining of this paper the subsampling in UV-CSXI and RV-CSXI is only performed using block-unblock coded apertures. The performance of CSXI, as well as that of the proposed reconstruction algorithm, is analyzed for each subsampling strategy and five different subsampling rates: 1 (full sampling), 4, 8, 16 and 32. Note the full sampling is achieved using the configurations shown in Fig. 2, while the subsampling rates greater than 1 are attained using the configurations shown in Fig. 3 by varying the number of blocking elements accordingly. The number of measurements acquired in each configuration is set to be the same such that the radiation dose is equivalent in all configurations for each subsampling rate. Table 2 shows the PSNR of each energy binned reconstruction for all subsampling strategies and subsampling rates. Note for RDB-CSXI, UVB-CSXI, and RVB-CSXI the PSNR values are the result of averaging 10 different random measurement selections in each case. It is observed that high-quality images at all energy bins can be attained using all the configurations and the proposed algorithm for the full-sampling scenario. It should be noted that the UDB-CSXI subsampling has the lowest performance compared to the other architectures for all subsampling rates. This is due to the fact that uniform sampling does not exploit the structure of the tomography sensing matrix, nor the structure of the filtered spectra in each case. Optimization of the coded apertures in the multiple settings is a topic of current work as it would lead to images with higher quality when the subsampling rate is low.

 

Fig. 9 (a) 3D-IRCADb thoracic-abdominal CT scan slice, (b) Segmented CT scan according to pixel material. (c-f) Original energy binned effective linear attenuation coefficient images. Reconstructions obtained with (g-j) UV-CSXI, and (k-n) RD-CSXI using one of 10 random masks at each view. The PSNR is shown for each reconstruction. (o-v) Normalized absolute error for each configuration.

Download Full Size | PDF

An additional phantom is generated based on the thoracic-abdominal CT scan, depicted in Fig. 9(a), which was obtained from the 3D-IRCADb database. To generate the multi-energy phantom, segments of the image are grouped according to the segmentation provided in [40]; as it can be seen in Fig. 9(b), the CT scan has the following materials: (1) lung, (2) heart (blood), (3) aorta (blood), (4) bone, (5) stomach (soft tissue), (6) adipose tissue, (7) air. The mass density for each pixel is estimated using a linear mapping of the grayscale value in the CT scan such that the resulting average mass density in the lungs is 0.32 $g/c{m}^3$. Using the mass density and the mass attenuation coefficients of the materials in the phantom obtained from the NIST tables [35], the CSXI projection data were simulated as described for the FORBILD thorax phantom. Figures 9(c)-9(f) depict the effective linear attenuation coefficients of the object for each energy bin. Note that an additional segment was added to the phantom representing iodine+blood in the heart area. Figures 9(g)-9(j) depict the reconstructions obtained for UV-CSXI (full sampling). In addition, a sampling scheme that uses only 10 different block/unblock masks is evaluated since this offers potential simplification for the RD-CSXI implementation. In this approach, the mask is still changed for every view, but there are fewer masks than views, so the masks are reused periodically (10 masks instead of 512). Figures 9(k)-9(n) depict the reconstructions for the proposed simplified architecture at all energy bins and the respective PSNR. Interestingly, note that the reconstruction quality is not degraded and the architecture is significantly simplified as evidenced by the PSNR and the normalized absolute error images (Figs. 9(o)-9(v)). In addition, note that in both the Forbild thorax phantom and the 3D-IRCADb scan the PSNR of the first energy bin is lower compared to the other energy bins particularly for low subsampling ratios, this is likely caused by the high fluctuation of the linear attenuation coefficients of the objects in the phantom in the first energy bin compared to other energy bins. This effect could be ameliorated by adding more filters with a K-edge in the first energy bin interval.

 

Fig. 10 (a) Effective linear attenuation coefficient images at energy bin k = 1. Reconstructions obtained using (b) UV-CSXI with $M=P=512$ and a single scan, and (c) the system proposed in [12] with $M=512$, $P=103$ and 5 different scans. The PSNR is shown for each reconstruction as well as zoomed versions of the highlighted areas in the phantom. The normalized absolute error is shown in each case.

Download Full Size | PDF

Table 3. Thorax DICOM reconstruction PSNR for UV-CSXI downsampling x8 vs K-edge filter system in [12] with $P=13$

View Table | View all tables in this article

In addition, to evaluate the performance of the CSXI system and the proposed reconstruction algorithm, reconstructions using the system developed by Rakvongthai et al. in [12] are performed and compared with the proposed method. In order to compare both systems, the radiation dose is set to be equivalent in both systems. For $M=P=512$ in CSXI, the number of angles in the system in [12] is set to $P=103$, Least Squares (LS) is used to obtain the energy binned sinograms and the energy binned attenuation coefficients are independently reconstructed using the 2D-FBP algorithm available in MATLAB (ifanbeam function) as described in [12]. The reconstructions of both systems are compared to the effective linear attenuation coefficients calculated as a weighted sum of the linear attenuation values obtained from the NIST database. The reference image at the first energy bin and the reconstructions of UV-CSXI and the system in [12] are presented in Figs. 10(a)-10(c), respectively. It is observed that high-quality images can be obtained using the CSXI system. On the other hand, the results obtained using the system proposed in [12] present numerous artifacts since FBP is not a suitable algorithm to solve the inverse problem when the number of view angles is limited. The performance of the limited view angle version of [12] would be improved by using an algebraic reconstruction approach with TV regularization. Table 3 includes the PSNR of the reconstructions for $P=13$ of a Thorax-DICOM phantom obtained using this latter approach, as well as FBP. Additionally, these reconstructions are compared with the CSXI system of an equivalent those, that is, the UVB-CSXI system with a downsampling rate of 8. It should be noted, however, that using TV regularization in the approach in [12] would not ameliorate the requirement of having to align F different consecutive scans for the system which in practice is often difficult.

5. Conclusion and future work

A spectral computed tomography system referred to as compressive spectral X-ray imaging (CSXI) which uses conventional X-ray imaging systems combined with K-edge coded aperture filters is proposed. CSXI has the potential to provide a cost-effective alternative to photon counting spectral CT techniques and a faster scanning solution with reduced X-ray radiation dose compared to traditional spectral CT imaging systems. In addition, as the system uses conventional integrating detectors, the spatial resolution that can be achieved in CSXI is higher compared with photon counting detector systems in general. Multiple CSXI subsampling architectures capable of reconstructing high-quality spectral datacubes using a single scan were analyzed. It was shown that the spectral datacube can be recovered from fewer measurements than the full-sampling scenario for all the proposed configurations. Furthermore, simulations showed an improvement of up to 7 dB in the reconstructions using the CSXI system compared to the system proposed in [12] when using FBP. The simulations were performed for a fan-beam geometry and the object was represented using a Kronecker product between a 2D Haar wavelet basis and a DCT basis, but CSXI is far more general and can be applied to other geometries and other sparse representations of the object as well. A three-phase reconstruction algorithm based on the estimation of the energy binned sinogram for each energy bin was proposed to reconstruct the energy binned linear attenuation coefficients. Other algorithms that further exploit the structure of the filtered sinograms could be used. Fabrication of the coded apertures and their integration with the K-edge filters for circular tomography is under development. Calibration procedures will be used in order to mitigate the mismatching errors that might occur. Furthermore, the angular collimation produced by the implemented coded apertures will be accounted for in the system matrix H. Based on the structure of the resulting sensing matrix, the optimization of the coded apertures of the CSXI configurations is a topic of current research.