Spectral computed tomography (CT) relies on the spectral dependence of X-ray attenuation coefficients to separate projection measurements into more than two energy bins. Such data can be used to unveil tomographic material characterization — key in national security and medical imaging. This paper explores a radical departure from conventional methods used in spectral imaging. It relies on K-edge coded apertures to create spatially and spectrally coded, lower-dose, X-ray bundles that interrogate specific voxels of the object. The new approach referred to as compressive spectral X-ray imaging (CSXI) uses low-cost standard X-ray integrating detectors and acquires compressive measurements, which enable the reconstruction of energy binned images from fewer measurements. Various spectral and spatial coding strategies for structured illumination are explored. Subsampling in CSXI is accomplished by either view angle spectral subsampling, spatial subsampling enabled by block-unblock coded apertures placed at the source or detector side, or both. The careful design of subsampling strategies, spectral filters, coded apertures, and their placement, are shown to be critical for the quality of tomographic image reconstruction. The forward imaging model of CSXI, which is a non-linear ill-posed problem, is analyzed and a multi-stage algorithm is developed to address the estimation of the energy binned sinograms from the integrating detector measurements. Then, an Alternating Direction Method of Multipliers (ADMM) is used to solve a joint sparse and low-rank optimization problem for reconstruction that exploits the structure of the spectral X-ray data cube.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
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 . 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 , or a set of sandwiched detectors capable of discriminating between low and high X-ray energies . 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 . 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 . 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  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  based on the Ross filter method . 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.  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 , spectroscopy  and more recently in X-ray CT monochromatization in medical imaging .
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.
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 . 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 . 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  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
When a polychromatic X-ray beam passes through an object, the X-ray attenuation at each energy, E, is obtained as , where is the linear attenuation coefficient of the object at position , is the X-ray source spectrum, is the energy-dependent detector response, and is the effective energy flux. For photon counting detectors and for conventional integrating detectors . 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 . Let 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 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 . The integrated gray-scale measurements in Eq. (1) cannot be used to reconstruct directly the energy binned images , as they do not provide spectrally resolved data . 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 . 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 . 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 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 . 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  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 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 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 sp, for , 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 , where is the linear attenuation coefficient of the filter f at the energy E, is the X-ray source intensity at energy E, and is the length of the intersection of the X-ray beam with the filter f. The latter is given by , where ρf is the thickness of the filter and ψj is the angle between the normal of the filter and the X-ray beam. Thus, the filtered measurement at the integrating detector element is given by:Eq. (3) cannot be rewritten as a set of linear equations in terms of the energy binned linear attenuation coefficients since they are a weighted sum of exponential functions of the spectral datacube. Consequently, the reconstruction problem is formulated as follows:
This minimization problem is non-linear and the set of measurements 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 . Note that the set of filtered projections measured by the integrating detector, in essence, builds up a set of F tiled subsampled sinograms , 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 , defined as , for (Fig. 4(b)). Unfortunately, and 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 which yields the effective linear attenuation coefficients using filter f, (Fig. 4(c)), and after the inverse log-normalization operation the estimate is obtained (Fig. 4(d)). These procedures are detailed in Section 3.1. The next phase, once the filtered intensities 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 . These steps will be detailed in Section 3.2. Finally, the effective linear attenuation coefficients at each energy bin 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.
3.1. Filtered sinograms inpainting
Each measurement 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 , which are then normalized to obtain the log-normalized sinograms as shown in Fig. 4(b). Then, each subsampled set of measurements is used to reconstruct the effective linear attenuation coefficients using regularization constraints to account for the missing data as depicted in Fig. 4(c). Subsequently, the effective linear attenuation coefficients can be re-projected to the sinogram space using the information of the CT system matrix H to obtain the missing points in 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 : , where I0 is the intensity of the X-ray source at the effective energy, I is the measured intensity on the detector, and 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:25]. As shown in , the under-determined system of equations can be solved as
While many algorithms exist to solve Eq. (6), the Two-Step Iterative Shrinkage/Thresholding (TwIST) algorithm  is used here to reconstruct and the Total Variation (TV) norm of is used as the regularization term. Then, the inpainted sinogram is attained using the complete system matrix to estimate the missing points, i.e. . The estimated intensity at the detector element j using the filter is then obtained as . Note, this process might change the values of at the locations that had been originally measured. Therefore, the original measurements are substituted into the estimated intensity sinograms; mathematically, , where Ωf corresponds to the set of values of j where the filter f was used.
3.2. Energy binned sinograms estimation
Given the estimated sinograms (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 energy bins , where Ef for is the K-edge energy of the filter f, , and ET 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 (E1) for Ce, 53.8 keV (E2) for Dy, and 57.5 keV (EF) 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 and is almost identical as shown in Figs. 6(b) and 6(c); and for the filter pair Er-Dy, the filtered spectra between and 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 and is almost identical. Generalizing for a set of F K-edge filters with balanced thicknesses, the energy spectra in the intervals and 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 , where , and is the monochromatic intensity at the detector element j at energy bin k, (see Appendix A). Given that the energy spectrum in the intervals and are very similar in value for all the filters, the recovery of these bands (k = 0 and ) is not possible. Let be the estimated filtered intensities, and be the energy binned intensities at the detector element. Then, the filtered intensities at the detector can be expressed as , where 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 detector, the remaining filtered intensities are estimates obtained by the process described in Section 3.1 and may contain errors. Thus, the pseudo-inverse of the matrix cannot be used directly to obtain the monochromatic intensities. In order to account for the noisy estimates, , a joint estimation of the monochromatic X-ray intensity sinograms using a gradient descent algorithm is used to solve the inverse problem:Fig. 5(b). Normalizing and taking the natural logarithm, the log-transformed sinograms for are obtained (Fig. 5(c)). Note that the estimated intensities for k = 0 and are not taken into account given that the recovery of these bands is not possible from the filtered measurements as previously explained.
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 at energy bin k, for . Using the monochromatic X-ray Beer-Lambert law, the energy binned intensities at the detector can be written as , where is the quasi-monochromatic intensity at energy k. Normalizing and discretizing the line integrals and the linear attenuation coefficients , the energy binned measurements (Fig. 5(c)) are given by the following linear model: , where is the vectorized linear attenuation coefficient of the object at the 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 . Here, a joint sparse and low-rank optimization method is used, but other approaches such as , 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 norm of the sparse representation of the datacube, and its nuclear norm. Furthermore, the tensor modeling previously proposed by Gao et al. in  and  is adapted to formulate the reconstruction problem as follows:
Vectorization: Let , be the operator that stacks the entries of a tensor in reverse lexicographical order into a -long column vector .
Unfolding: For a tensor , the matrix unfolding denoted by is defined by the mapping of: . The resulting matrix is a 2D matrix with Qn rows and columns .
Folding: The operator is the inverse operator of unfolding .
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 , 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 is a multi-dimensional matrix, where the first two dimensions are spatial, the third is the energy dimension, and . The spectral CT model can thus be generalized as:Eq. (8). As previously mentioned, regularization for denoising is used to reconstruct the datacube. Let be represented by , where is the sparse tensor representation of the object, and is the representation basis. Then, the sparse representation of the object can be reconstructed by solving the optimization problem .
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 [10, 18] as:31] is adapted to solve the inverse problem as it provides the necessary tools to split the optimization problem into small convex optimization problems that can be solved using simpler algorithms. When the augmented Lagrangian is calculated, two intermediate variables to separate the , and nuclear norm components are introduced. Thus, Eq. (9) is transformed as:
The problem is then split into three different sub-problems summarized in Algorithm 1. Step 2 of the algorithm is an minimization problem which is solved using the conjugate gradient (CG) algorithm with a fixed step size, and as introduced in  the solution for Step 3 and Step 4 are respectively given by:10], and are the matrices obtained by applying the unfolding operator to the tensors and , respectively and is the shrinkage operator with parameter δ operating in each component of as:
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 . Hence, in order to find a suitable initialization point for the ADMM algorithm, the energy binned images are reconstructed from the corresponding quasi-monochromatic sinograms by solving the following minimization problem independently for each energy bin:Eq. (14) is solved independently for each energy bin using the TwIST algorithm  and the resultis used as the initial value, , for the joint reconstruction algorithm detailed in Algorithm 1.
4. Simulation results
ACT configuration with a flat 1D detector strip with elements, a fan beam X-ray source, and view angles was used in the simulations. The scanned object was a 256 × 256 modified Forbild thorax phantom  generated using the CONRAD software . 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 . 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 cm with , to match the full sampling scenario described by Jorgensen et al in . 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 , 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 respectively, as detailed in Table 1. The X-ray filtered spectra were simulated at 80 keV using the Spektr software  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.
The projection data were simulated using the ASTRA tomography toolbox  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 . 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.
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 , where is the Kronecker product, 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 . The peak signal to noise ratio (PSNR) is used to evaluate the reconstructions since it is suitable for comparing restoration results . 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.
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 ; 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 . Using the mass density and the mass attenuation coefficients of the materials in the phantom obtained from the NIST tables , 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.
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  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 in CSXI, the number of angles in the system in  is set to , 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 . 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  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  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  would be improved by using an algebraic reconstruction approach with TV regularization. Table 3 includes the PSNR of the reconstructions for 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  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  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.
Appendix A : Transmission matrix B weights
The integrating detector measurements in Eq. (2) can be rewritten as , where corresponds to the X-ray intensity that would be measured by an energy selective detector in the energy bin, using filter f and it is given by:13]; thus, the linear attenuation coefficient in the energy bin can be approximated by the average attenuation coefficient in the energy bin and Eq. (15) can be simplified as:
National Science Foundation (NSF) (CIF 1717578); University of Delaware UNIDEL Grant.
1. B. J. Heismann, B. T. Schmidt, and T. G. Flohr, Spectral Computed Tomography(SPIE Press Book, 2012). [CrossRef]
2. B. Krauss, B. Schmidt, and T. G. Flohr, Dual Source CT(Springer, 2011), pp. 11–20.
3. R. Carmi, G. Naveh, and A. Altman, “Material separation with dual-layer CT,” in IEEE Nuclear Science Symp. Conf. Record, 2005, vol. 4 (2005), pp. 3 pp.–1878.
4. J. P. Schlomka, E. Roessl, R. Dorscheid, S. Dill, G. Martens, T. Istel, C. Baumer, C. Herrmann, R. Steadman, G. Zeitler, A. Livne, and R. Proksa, “Experimental feasibility of multi-energy photon-counting K-edge imaging in pre-clinical computed tomography,” Phys. Medicine Biol. 53, 4031 (2008). [CrossRef]
5. A. P. Cuadros and G. R. Arce, “Coded aperture compressive X-ray spectral CT,” in 2017 Int. Conf. on Sampling Theory and Applications (SampTA), (2017), pp. 548–551.
6. I. Glass, A. P. H. Butler, P. H. Butler, P. J. Bones, and S. J. Weddell, “Physiological gating of the MARS spectral micro CT scanner,” in 28th Int. Conf. on Image and Vision Computing New Zealand (IVCNZ), (2013), pp. 317–322.
7. C. T. Badea, D. P. Clark, M. Holbrook, M. Srivastava, Y. Mowery, and K. B. Ghaghada, “Functional imaging of tumor vasculature using iodine and gadolinium-based nanoparticle contrast agents: a comparison of spectral micro-ct using energy integrating and photon counting detectors,” Phys. Medicine Biol. (2019). (posted February 1st 2019, in press). [CrossRef]
9. K. Taguchi, M. Zhang, E. C. Frey, X. Wang, J. S. Iwanczyk, E. Nygard, N. E. Hartsough, B. M. W. Tsui, and W. C. Barber, “Modeling the performance of a photon counting X-ray detector for CT: Energy response and pulse pileup effects,” Med. Phys. 38, 1089–1102 (2011). [CrossRef] [PubMed]
10. L. Li, Z. Chen, W. Cong, and G. Wang, “Spectral CT modeling and reconstruction with hybrid detectors in dynamic-threshold-based counting and integrating modes,” IEEE Transactions on Med. Imaging 34, 716–728 (2015). [CrossRef]
11. Hamamatsu Handbook, “X-ray detectors-chapter 09,” https://www.hamamatsu.com/resources/pdf/ssd/e09_handbook_xray_detectors.pdf. Accessed March 7, 2019.
12. Y. Rakvongthai, W. Worstell, G. E. Fakhri, J. Bian, A. Lorsakul, and J. Ouyang, “Spectral CT using multiple balanced K-edge filters,” IEEE Transactions on Med. Imaging 34, 740–747 (2015). [CrossRef]
13. P. A. Ross, “Polarization of X-rays from a tungsten target,” J. Opt. Soc. Am. 16, 375–379 (1928). [CrossRef]
14. W. Bol, “The use of balanced filters in X-ray diffraction,” J. Sci. Instruments 44, 736 (1967). [CrossRef]
15. W. Miao, Y. Ding, K. Sun, D. Lai, Z. Chen, H. Wang, J. Cheng, Z. Zheng, and H. Peng, “Design of the Ross pairs for soft X-ray spectrometer,” in IEEE Int. Conf. on Plasma Science (Cat. No.02CH37340), (2002), p. 161.
16. M. Saito, “Quasimonochromatic X-ray computed tomography by the balanced filter method using a conventional X-ray source,” Med. Phys. 31, 3436–3443 (2004). [CrossRef]
18. H. Gao, H. Yu, S. Osher, and G. Wang, “Multi-energy CT based on a prior rank, intensity and sparsity model (PRISM),” Inverse Probl. 27, 115012 (2011). [CrossRef]
19. J. Prince and J. Links, Medical imaging signals and systems (Pearson Prentice Hall, 2006).
20. J. Chu, W. Cong, L. Li, and G. Wang, “Combination of current-integrating and photon-counting detector modules for spectral CT,” Phys. Medicine Biol. 58, 7009 (2013). [CrossRef]
21. B. D. Arhatari and T. E. Gureyev, “Elemental contrast X-ray tomography using Ross filter pairs with a polychromatic laboratory source,” Sci. Reports 7, 218 (2017).
22. P. Ross, “Polarization of X-rays,” Phys. Rev 28, 425 (1926).
23. P. Kirkpatrick, “On the theory and use of Ross filters,” Rev. Sci. Instruments 10, 186–191 (1939). [CrossRef]
24. Y. Kaganovsky, D. Li, A. Holmgren, H. Jeon, K. P. MacCabe, D. G. Politte, J. A. O’Sullivan, L. Carin, and D. J. Brady, “Compressed sampling strategies for tomography,” J. Opt. Soc. Am. A 31, 1369–1394 (2014). [CrossRef]
25. E. J. Candes and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. 25, 21–30 (2008). [CrossRef]
26. J. M. Bioucas-Dias and M. A. T. Figueiredo, “A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration,” IEEE Transactions on Image Process. 16, 2992–3004 (2007). [CrossRef]
27. K. Kim, J. C. Ye, W. Worstell, J. Ouyang, Y. Rakvongthai, G. E. Fakhri, and Q. Li, “Sparse-view spectral CT reconstruction using spectral patch-based low-rank penalty,” IEEE Trans. on Medical Imaging 34, 748–760 (2015). [CrossRef]
28. L. Li, Z. Chen, G. Wang, J. Chu, and H. Gao, “A tensor PRISM algorithm for multi-energy CT reconstruction and comparative studies,” J. X-Ray Sci. Technol. 22, 147–163 (2014).
29. D. Kressner, “Low-rank tensor techniques for high-dimensional problems,” http://www2.mpi-magdeburg.mpg.de/mpcsc/events/trogir/slides/tensorlectures.pdf. Accessed March 1, 2018.
30. L. D. Lathauwer, B. D. Moor, and J. Vandewalle, “A multilinear singular value decomposition,” SIAM J. Matrix Analysis Appl. 21, 1253–1278 (2000). [CrossRef]
31. S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Found. Trends Mach. Learn. 3, 1–122 (2011). [CrossRef]
32. K. Hamalainen, L. Harhanen, A. Hauptmann, A. Kallonen, E. Niemi, and S. Siltanen, “Total variation regularization for large-scale X-ray tomography,” Int. J. Tomogr. Simul. 25, 1 – 25 (2014).
33. K. Sourbelle, “Forbild thorax phantom,” http://www.imp.uni-erlangen.de/phantoms/thorax/thorax.htm. Accessed March 1, 2018.
34. A. Maier, H. G. Hofmann, M. Berger, P. Fischer, C. Schwemmer, H. Wu, K. Muller, J. Hornegger, J. H. Choi, C. Riess, A. Keil, and R. Fahrig, “Conrad-a software framework for cone-beam imaging in radiology,” Med. Phys. 40, 111914 (2013). [CrossRef]
35. J. Hubbell and S. Seltzer, “Radiation and biomolecular physics division, MPL NIST,” https://physics.nist.gov/PhysRefData/XrayMassCoef/tab3.html. Accessed March 1, 2018.
36. J. S. Jorgensen, E. Y. Sidky, and X. Pan, “Quantifying admissible undersampling for sparsity-exploiting iterative image reconstruction in X-ray CT,” IEEE Transactions on Med. Imaging 32, 460–473 (2013). [CrossRef]
37. J. H. Siewerdsen, A. M. Waese, D. J. Moseley, S. Richard, and D. A. Jaffray, “Spektr: A computational tool for X-ray spectral analysis and imaging system optimization,” Med. Phys. 31, 3057–3067 (2004). [CrossRef] [PubMed]
38. W. Xu, F. Xu, M. Jones, B. Keszthelyi, J. Sedat, D. Agard, and K. Mueller, “High-performance iterative electron tomography reconstruction with long-object compensation using graphics processing units (gpus),” J. Struct. Biol. 171, 142 – 153 (2010). [CrossRef] [PubMed]
39. L. Galvis, H. Arguello, and G. R. Arce, “Coded aperture design in mismatched compressive spectral imaging,” Appl. Opt. 54, 9875–9882 (2015). [CrossRef]
40. L. Soler, “3d image reconstruction for comparison of algorithm database: A patient specific anatomical and medical image database,” http://www.ircad.fr/research/3d-ircadb-02/. Accessed Nov 1, 2018.